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SUBSONIC ANNULAR WING THEORY WITH APPLICATION 
TO FLOW ABOUT NACELLES 

By Michael J. Mann 
Langley Research Center 

SUMMARY 

Jet engine nacelles have generally been designed with the assumption that the intake 
and nozzle can be treated separately. However, advancements in engine development for 
subsonic cruise airplanes have tended toward high -bypass -ratio fan engines. This type of 
engine usually has a short fan nacelle which necessitates analyzing the complete nacelle. 
These nacelles can be treated as annular wings on which the circulation developed deter- 
mines both the internal and external flow. 

A method has recently been developed for calculating the flow over such a nacelle at 
zero angle of attack and at subsonic Mach numbers. The method makes use of annular 
wing theory and boundary -layer theory and has shown good agreement with both experimen- 
tal data and more complex theoretical solutions. The method permits variation of the 
mass flow by changing the size of a center body. 

INTRODUCTION 

Jet engines have generally been designed with the assumption that the inlet and nozzle 
can be designed separately. Advancements in engine technology have recently tended 
toward high -bypass -ratio fan engines with short fan nacelles. For these nacelles, the com- 
plete nacelle must be examined. A complete nacelle can be viewed as an annular wing, the 
thickness and camber of which determine the internal and external flow. 

Several methods have been developed for calculating the flow over nacelles or annular 
wings. The Douglas -Neumann method (ref. 1) has been used to predict accurately the pres- 
sures on the surface of a nacelle inlet. Trulin and Iversen (ref. 2), Geissler (ref. 3), and 
Young (ref. 4) have recently published methods which use surface singularities and show 
good agreement between theory and experiment for complete nacelles. Geissler*s method 
includes the effects of angle of attack. 

Young (ref. 5) also examined the accuracy of a method which distributes the singu- 
larities along a-mean cylinder of constant diameter. However, Young did not present the 
details of the theory. Agreement between the theory and experiment is good at low Mach 
numbers but is inadequate at the higher Mach numbers. 



Mascitti (ref. 6) has developed a rapid numerical technique for calculating the sub- 
sonic flow over planar and axisymmetric profiles at zero angle of attack. The method 
solves the exact incompressible potential -flow equations written in terms of the stream 
function. 

Belotserkovskii (ref. 7) and Weissinger (ref. 8) have developed much of the theory 
of annular wing flow. Their application of the theory has been the determination of over- 
all force and moment characteristics rather than a study of surface pressure distribution. 
However, much of their theoretical work has been useful in the present study. 

Hess (ref. 9) discusses some of the early solutions for axisymmetric flow which 
made use of vortex rings on the surface. Keith and others (ref. 10), and Grossman and 
Moretti (ref. 11) have developed numerical methods for solving the exact inviscid flow 
over two-dimensional and axisymmetric nacelles at subsonic and transonic speeds. Ref- 
erence 10 uses a stream-tube curvature relaxation technique, and reference 11 uses a 
time -dependent technique. However, the technique of reference 11 was only successful 
for subsonic flow. To the author's knowledge, only the present paper and references 3 
and 10 have included the effects of the boundary layer. 

The purpose of the present study is twofold. First, an inviscid theory, which uses 
vortex rings to represent the camber and a combination of source rings and vortex rings 
to represent thickness has been developed. Three-dimensional sources are distributed 
along the center line to control the mass flow through the wing. This method appears to 
have some advantage over each of the other methods in terms of either simplicity, a 
reduction of computer core storage or central processor time requirements, ease of exten- 
sion to more complex flows, accuracy, or quantities computed. Secondly, an analysis of 
the importance of including the displacement effect of the boundary layer in the calcula- 
tions has been made. A boundary -layer iteration procedure has been included in the com- 
puter program. 

Results of the present theory have been compared with experimental data obtained 
in tests on four nacelles and good agreement has been found over a range of mass -flow 
ratios and Mach numbers. 


SYMBOLS 

internal -flow cross-sectional area, m 


2 (ft2) 


wing chord, m (ft) 


pressure coefficient. 


P - P. 


2 



d 

D 




E 

gq(4) 

H 

K 


k 


k’ 


M 


OO 


m 


N 

P 

QU) 

q(4) 

r 

^00 

tU) 


intake diameter equals highli^t diameter less twice inlet lip radius, m (ft) 

average of minimum inner diameter and maximum outer diameter of wing; 
when subscripted, refers to diameter of wing, m (ft) 

unit vectors in meridional, r- and x-directions, respectively 

p ^/2 / 2 ^ 

complete elliptic integral of second kind, \ vl - k sin v du, dimensionless 

dimensionless correction term for (see eq. (32)) 
shape factor, b*/Q 

pV2 . 

complete elliptic integral of first kind, \ ■ — do. dimensionless 

*^0 I ,2 . 2 

- k sm V 

modulus of complete elliptic integrals (see eq. (11)), dimensionless 
complementary modulus, yl - k^ 
free -stream Mach number 

mass flow through annular wing, kg/ sec (slugs/ sec); iteration number or 
summation index 

number of chordwise divisions for numerical solution 
static pressure, N/m^ (lbf/ft2) 

thickness function defined by equation (29), dimensionless 
source strength per xinit axial distance, m/ sec (ft/ sec) 
radial coordinate attached to wing, m (ft) 
free-stream dynamic pressure, PoqV^/ 2, N/m^ (lbf/ft2) 
wing thickness, router - ^inner’ “ 
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Uq(^ - dimensionless kernel defined by equation (31) for = 1 

dimensionless kernel defined by equation (30) for Pq = 1 
V total velocity, free stream plus perturbation, m/ sec (ft/ sec) 


free -stream velocity, m/sec (ft/sec) 



perturbation velocity components in x- and r-directions, respectively, induced 
at m/sec (ft/sec) 


«(« - Vo) 
'’(« - «o'Po) 


dimensionless perturbation velocity influence functions for x- and 
r-directions, respectively; these influence functions are for velocity 
induced at ^q,Pq by a ring singularity at ^,p >= 1 or by a point 
singularity at ^,p = 0 


X inlet length, m (ft) 


X, y, z, r coordinates attached to wing (see fig. 1 and app. A for positive directions), 
m (ft) 


Q!q slope of mean camber surface with respect to x-axis (wing centerline), radians 

slope of symmetrical thickness profile with respect to x-axis, (dt/dx)/2 
radians 


Fq dimensionless bound vortex strength of camber solution, 7 ^ 

dimensionless boxmd vortex sheet strength of camber solution, y^yv^o^ 
y^ bound vortex sheet strength of camber solution, m/sec (ft/sec) 

y^ vortex sheet strength for thickness, m/ sec (ft/ sec) 

y ratio of specific heats 
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inner minus outer of ( ); A| = dimensionless panel length 


boundary-layer displacement thickness, m (ft) 


small positive number 


momentum thickness, m (ft) 


aspect ratio, D/c 


mass-flow ratio (see eq, (41) 


dimensionless distance along axis. 


dummy variable form of | 


dimensionless radius, — gas density, kg/m^ (slugs/ ft^) 

U/ Ci 

dummy variable 


digamma function (see app. B) 


Superscripts: 


incompressible quantity except and k' 


unit vector; velocity influence function 


vector, except gq(^) 

critical Mach number condition; boundary -layer displacement thickness 


Subscripts: 


dimensional bound vortex sheet strength of camber solution 


free stream 


center body 
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exp experimental 

h highlight station (most forward point on inlet lip, x = c/2 in fig. 1) 

i location of vortex and source singularities 

j location of point at which velocity is induced 

LE leading edge 

m iteration number (see eq. (44)) 

max maximum 

min minimum 

o point at which velocity is induced in effect due to camber only without 

angle -of -attack or pitching-motion effects in a^, and r^; zeroth 
order in Fourier expansions in also see equation (44) 

q due to source ring of thickness solution 

r radial direction 

TE trailing edge 

t thickness solution; tangential 

th theory 

w wing surface 

X wing axial direction 

y due to vortex ring 

axial location 

pQ radial location 
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THEORETICAL SOLUTION 


The subsonic flow over an axisymmetric nacelle or annular wing with a center body 
of revolution will be determined by use of small perturbation theory. This problem is 
analogous in many ways to airfoil theory and biplane theory. Linearization of the problem 
permits computation of the compressibility effects by transforming the wing coordinates 
X and r into an equivalent incompressible wing x’ and r' by using the Gbthert 
transformation 


X* = X 

(1) 

r’ = |3r 

(2) 


where /3 = yl - M^. The incompressible flow over this equivalent wing is then computed, 
and the compressible perturbation velocities u and v in the x- and r-directions, 
respectively, are found by 


u = u’/i3^ 

(3) 

<J 

11 

< 

(4) 


while the free-stream velocity remains unchanged. Compressible flow equations 
can then be used to determine the pressure distribution. (All primes are dropped from 
coordinates and velocities from this point on, and these variables are understood to be 
incompressible, unless otherwise stated.) 

Also as a consequence of the linearization of the equations of motion and the boundary 
conditions, the problem can be broken up into wing camber, wing thickness, and center-body 
portions. Furthermore, each of these problems can be solved by superimposing elementary 
singular solutions of Laplace’s equation, such as vortexes and sources. The linearized 
boundary condition on the wing surface at the point ?q»Pq is 



where dr/dx is the local slope of the wing and and are dimensionless x- and 
r-coordinates. The minus sign arises because the positive x-direction has been chosen 
to be upstream (fig. 1). This boundary condition is written in two parts, one part for the 
camber and center body- and one part for the thickness 
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(6) 



V 


a 




(7) 


where and are the slopes of the mean camber line and the symmetrical thick- 
ness, respectively. The subscripts y, b, q, and t,y refer to camber vorticity, center 
body, thickness source, and thickness vorticity, respectively. Two terms arise in the 
thickness equation because sources alone will not produce a symmetrically thick wing. 

This point is discussed further when the thickness solution is outlined. When plus and 
minus signs are used, the upper and lower signs refer to the outer and inner wing sur- 
faces, respectively. 

For design purposes it is convenient to separate the problem into its camber and 
thickness components. Thus, all singularities (except the center body) are placed on a 
cylinder of constant diameter D as shown in figure 1 and the boundary conditions are 
satisfied at control points on this same cylinder (coxmterpart of the planar wing approxi- 
mation). A correction for approximating the surface velocities by velocities on this cylin- 
der, which is called the Riegels factor, is discussed subsequently. 

The strengths of the singularities are foimd by solving the boundary conditions given 
by equations (6) and (7). The solution of these boundary conditions involves evaluating 
improper integrals. The technique used is simply to replace the integrals by summa- 
tions, thereby replacing the continuous vortex and source sheets and the source line by 
distributions of singularities. In the camber solution, proper placement of the vortex 
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Figure 1.- Coordinates, geometry, and singularity and control-point location. 



singularities and control points will satisfy both the principal value of the integral and the 
Kutta condition (ref. 7, p. 106). The same lattice arrangement is used for the sources as 
for the vortexes. 

Once the velocity expressions are substituted into equations (6) and (7), a set of 
linear equations results with xmknown singularity strengths. The vortex strengths are 
found by solving these equations by the method of successive orthogonalization. The 
source strengths for the wing thickness and the center body are proportional to the local 
body slopes. 

In the following sections the major equations are given. The camber— center-body 
solution, the thickness solution, the complete wing solution, and the boundary -layer calcu- 
lation are discussed. 


Camber— Center-Body Solution 

The flow to be calculated is subsonic flow over an axisymmetric mean camber line 
at zero angle of attack. Along the center line of this infinitesimally thin wing is mounted 
a body of revolution. The mass flow inside the camber surface is varied by varying the 
size of the center body. 

In order to obtain a unique solution for the flow over a camber surface, the basic 
flow conservation equations must be satisfied, boundary conditions on the body and at 
infinity must be satisfied, and the circulation specified (accomplished by the Kutta condi- 
tion). As mentioned previpusly, these conditions can be met by distributing a continuous 
vortex sheet along a cylinder the length of which is equal to the wing chord. The influence 
of the center body can be found by use of a source line along the wing axis. 

In this study the influence of the wing on the center body is ignored. There are 
several reasons for ignoring this interference. First, owing to the wide variation in inter- 
nal structure, such as stings, struts, rakes, and so forth, no attempt was made to simulate 
any particular internal structural geometry. Secondly, inclusion of the mutual interference 
would greatly increase the storage requirements of the computer program (see app. A), 
thereby losing one of the attractive features of this method. Furthermore, the primary 
objective has been to obtain an accurate solution for the external flow. Thus, a simple 
body of revolution on the wing axis was selected to represent the internal parts. The size 
of this body is selected to either match the desired mass flow or match the internal block- 
age (and compute the mass flow). Both methods have shown good success as is demon- 
strated in the next section. 

The vortex sheet strength is found by first substituting the expressions for the vor- 
tex velocities (refs. 7 and 12) and source velocities (refs. 12 and 13) into the boundary 
condition, equation (6). When this is done, the result is 
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where the influence functions are defined by 



and K and E are complete elliptic integrals of the first and second kinds with modulus 





The variables ^ and p are the dimensionless x- and r-locations of singularities 
(p = 1), and ^q,Pq is the point at which the velocity is calculated (^Pq = l) • The vari- 
able X is the wing aspect ratio; ygCl) is the vortex strength and is related to the usual 
definition of vortex sheet strength by 

y+ = ^co%ro 

where is proportional to Ap/pV^Q. 

The subscript o on y^ and is used in reference 7 to indicate camber effects 
without angle of attack or pitching motion. For the sake of consistency this notation is 
retained. When angle -of-attack effects are included, the vortex sheet is no longer of con- 
stant strength around its circumference and a trailing vortex system must also be included 
(ref. 7). The camber effect alone will give a nonzero radial force; however, owing to 
symmetry the net lift is zero. 

Equation (8) is solved numerically by replacing the continuous vortex sheet by N 
equally spaced vortex rings. The source line along the center-body axis is replaced by 
three-dimensional point sources at the same axial locations as the vortex rings. ^ The 

^Center-body length can be less than the wing chord c. The case of a semi-infinite 
center body can be simulated by a constant body radius to downstream infinity because the 
local source strength is zero where dp^jd^ = 0. 
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singularity location is i and the control-point location is j. When the quantity 
is replaced by 


I'd!' 


= Voi 


equation (8) becomes 




a .r -V 
01 01 y 


(«i - = 



( 12 ) 


Because the vortex rings and control points are located on a constant -diameter right 

p 

circular cylinder of diameter D, then = p^j = 1 in equation (12), Each circular 
panel has a chord length of c/N, The vortex ring is located at one -fourth of this length 
from the panel leading edge and the control point is located at three -fourths of this length 
from the panel leading edge (fig. 1). 


As previously mentioned, the influence of the wing on the body is neglected. Thus 
the body source strengths are given by the body slopes and everything on the right-hand 
side of equation (12) is known. Equation (12) is then a set of N simultaneous linear equa- 
tions with N imknown vortex strengths ^see definition of CAM in app. A for dis- 

cussion of case when = 6j. 

The set of equations resulting from equation (12) was solved for the unknown vortex 
strengths F^. by the successive orthogonalization procedure outlined in reference 1. 

This procedure gave the same results as an alternate matrix inversion technique but has 
a great advantage on computer storage requirements. The method of reference 1 requires 
storage for a single N-dimensional array, the elements of v^ich are ~ ^oj’^)*^oi 

as opposed to other matrix techniques which require N- by N-array storage. 


Once the vortex and source strengths are computed, the perturbation velocities can 
be found at each location on the cylinder p^. = 1 by use of the following expressions 



V O' .F . 

V 01 01 

Voo 

^ 2ir 


i=l 


Q! .F . 
V 01 01 , 

VoO 

2t7 

i = 1 


a fl- - I .,l'l ± 


(13) 


(14) 


2 

The radius p of a singularity ring always has a value of 1 in this paper; however, 
the radius of the point at which the velocity is computed p varies between 0 and 1 in the 
mass -flow computation but has a value of 1 elsewhere. 
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( 15 ) 



( 16 ) 


where the influence functions are 



(17) 

(18) 

(19) 


( 20 ) 


3 

For points off the cylinder = 1 (necessary when computing mass flow), equa- 
tions (13) to (16) apply with the ±term dropped in equation (13) and using the influence 
functions 



^The variables ^ and p locate a singularity, and is the point at which 

the velocity is calculated. When going from integrals and continuous distributions to sum- 
mations and discrete distributions, it is necessary to indicate these quantities by ^.,Pj 

and ^Qj,PQj, respectively. When p^. = 1, the location loj’^oj ^ ^ control point; but 

when Pqj f 1, then 1® meant to indicate any arbitrary point. 
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^y(^i ■ ^oj’^oj ) ^ equations (9) and (10), Equation (13) 

is derived in appendix B by using the influence function in equation (21). 

A linear interpolation of the F^. was used to find the F^j. A check of the F^. 
distributions showed that the distribution approaches zero at the trailing edge, thereby 
satisfying the Kutta condition. Note that the pressure distribution over the camber sur- 
face depends on three quantities; camber surface slopes, aspect ratio, and mass flow 
(controlled by the size of the center body). 

Thickness Solution 

The flow to be calculated in this case is the subsonic flow over an annular wing with 
a symmetrical airfoil section at zero angle of attack. The flow must satisfy the same 
basic equations and conditions as the camber solution. The tangent flow condition is now 
equation (7). 

In the case of planar problems, thickness can be accounted for by a distribution of 
sources on a plane. Such a distribution has the necessary symmetry property and discon- 
tinuous change of sign in the vertical velocity for the planar wing thickness. However, if 
two planar source sheets are placed in close proximity to form a biplane, the interference 
between the two sets of source distributions will cause a mutual induction of camber. The 
vertical velocity at some point on either sheet will no longer depend solely on the local 
source strength but will also depend on the complete source distribution of the other sheet. 
This interference problem also exists for an annular wing. The induced camber effects 
can be canceled by a vortex distribution, the strength of which is adjusted to just cancel 
the interference effects on the wing. Then the total source and vortex radial velocity will 
have the necessary discontinuous property across the singularity sheet. 

Note that, because of the interference effects, an annular wing with thickness only at 
zero angle of attack has a negative radial force; that is, it tends to collapse. The internal 
flow area reduction due to thickness accelerates the flow much more than the same flow 
area reduction externally. Figure 2 shows the pressure distributions for an annular wing 
with thickness only and the same wing with enough camber to make the inside diameter 
constant. The radial force changes sign when geometric camber is added, and with the 
proper amount of geometric camber the inner and outer pressures can be equalized. 

Because the annular wing with thickness only has the effect of camber in its flow, 
namely a nonzero radial force, it would be convenient to refer to this effect as an aero- 
dynamic camber. Then the variations of the mean camber line from the chord could be 
referred to as geometric camber. 

The radial velocity induced at and Pq = 1 by continuous distribution of 
source rings of diameter D (radius p = 1) over an axial length of c is (see app. B) 
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(a) 4.3 -percent thickness and zero camber; /J. = 0.87. 

Figure 2. - Pressure distribution on an annular wing without and with 

camber at = 0.50. 




!al!2d]= - 1 ,i)di 

* 2V«, J 1 4irV-, 


where q(^) is the source strength at ^ and the influence function - ^^,1 j is 

related to the vortex influence function for the velocity u by 

= (24) 

The integral term in equation (23) is the induced camber due to interference effects. 
This integral can be canceled by a vortex distribution the strength of which is determined 
by 

1 


2J1 


= JliiLe - 4 l)d^ 
J 1 AirV^ q'^ o’ / 


where y^.(^) has the dimensions of y^, and ^ is computed by using equation (18). 
When this is done, equation (23) can be written as 

.4(0 


The flow constructed in this way will satisfy the required boundary condition in equa- 
tion (7). Thus from equations (7) and (26) the source strength is determined by 

l(5o)= -2V»“t(«o) (22) 

where the boundary condition is being satisfied at points ^ Iq’^o “ ^)' 

The solution of equation (25) for the unknown vortex strength is the same type of 
numerical problem found in the camber solution. An alternate method of finding ) 
has been outlined by Weissinger (ref. 8). By examination of the boundary conditions for 
a vortex flow, the vortex strength necessary to cancel the integral of equation (23) is 
determined by 
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X 


i^(S) 


^o) 


( 28 ) 


where 


Q({) 


VcX5 ^ 

(29) 


X 





^°|[ 2 (l + k2)E - (2 + k2)^ 

2 ll-tolL - 

(30) 

Uo(« 



(31) 


and t($) is the wing thickness and 


iq(^) 


rt(^) 


9.^ . 9 


is related to the vortex strength by 


(32) 


( The symbol ~Y^(i)/yoo is called gq(0 in ref. 8. ) The kernel functions and 
are continuous and antisymmetrical. The solutions of equations (25) and (28) for yj.(^) 
should require about the same computing time; equation (28) was chosen for this study. In 
numerical form, equations (28) and (29) yield 

N ^ N 

i= 1 i= 1 


«o(«l - «ofl) 


^i ^oj 


X 


(33) 


The source and vortex-ring locations ^ |j,p. = 1 ) and control-point locations 
( P .,p . = 1 ^ are the same as in the camber solution. The same method of solution of 
linear equations discussed previously was used for equation (33). 

Once the source and vortex strengths are computed, the perturbation velocities can 
be computed at any point (^oj’^oj) Equations (13) and (14) give the vortex velocities u^ ^ 
and v^^ when is replaced by (yn/^oo)^^ and similarly for Recall 

that for points off the cylinder p^. = 1, the ± term in equation (13) must be dropped. The 
source velocities at points on the cylinder p^. = 1 are computed from (refs. 8 and 12) 
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( 34 ) 





where the influence function ^ is given by equation (24) and u 


(35) 


"q(«l - C 

^O] 

Equation (35) can be obtained from equations (23) and (27). For points off the cylinder 
defined by = 1 (necessary when computing the mass flow), equations (34) and (35) 
apply with the term in equation (35) dropped and by using the influence fimctions 




Note that the thickness solution depends on two quantities : thickness distribution 
t(^) and aspect ratio. 

Pressure Distribution and Mass Flow 

Once the singularity strengths are determined, the velocities due to camber, center 
body, and thickness can be added to give the total surface velocity and pressure distribu- 
tion. The internal flow velocities can be integrated over a surface perpendicular to the 
center line at some ^-location to calculate the mass flow through the annular wing. 
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In regions of high body curvature, such as the nose, approximation of the surface 
velocity by the velocity on a constant -diameter cylinder is not satisfactory. A correction 
for this approximation called the Riegels factor, is shown by Weber (ref. 14) to improve 
the calculated velocities in the nose region of a two-dimensional airfoil. 

The Riegels factor can be developed from geometry as follows. By writing the 

4 

equation for the wing surface as 
f(x,r) = r - r^(x) = 0 

the unit tangential vector in the meridional direction is 



where and are unit vectors in the x- and r-directions, respectively. Writing 
the total compressible surface velocity as 



the total tangential surface velocity is 



Consistent with linear theory the total surface velocity is approximated by 



The quantity in the denominator is the Riegels factor. Equation (39) shows that the 
total velocity on the wing surface is foimd by computing the total x-velocity on the cylin- 
der p = 1 and multiplying this result by the Riegels factor. 

^All velocities from now on are compressible and coordinates are for the actual 
(untransformed) geometric shape. 
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The surface pressure distribution is computed by making compressibility and Riegels 
factor corrections to the incompressible velocities and calculating the pressure from 

r 


= 


JL. 


yM„ 




1 + 


2 ^ r - 1 


V. 


(40) 


j 


where y is the ratio of specific heats. 


The mass flow is computed in the form of a ratio 


4 = 


m 

Poo'^lj^oO 



(41) 


where m is the mass flow throu^ the annular wing, V t is the total x-velocity at some 

^o 

location, and is a reference area based on the highlight diameter (diameter at 
X = c/2 in fig. 1). The area is the internal cross-sectional area of the wing at 

less the local center-body cross-sectional area. The velocity ratio is computed by divid- 
ing the flow area into small annuli and calculating the mass flow through each annulus and 
summing the result as follows (assuming density constant over the cross section) 


_ V ^(^o’^oj) ^(^o’^oj) 
Vco ■ f A^^ 


(42) 


The V ^^o’^o j ) total x-velocities computed at a fixed and the particular . 

location. The annulus area is at the IqjPqj location. The density ratio is 

calculated from the one -dimensional relation 



y - 1 


(43) 


Boundary -Layer Solution 

The boundary layer was approximated by a turbulent boundary layer over the entire 
chord length. The incompressible turbulent boimdary-layer method of Truckenbrodt 
(ref. 15) was used. John B. Peterson, Jr., of NASA Langley Research Center developed 
the necessary equations for computing the boimdary layer and programed these equations 
for the boundary -layer subroutine used in this study. 
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The inviscid velocity distribution at the trailing edge is greatly in error since a 
stagnation point is computed there. Hence, the velocity was extrapolated to the trailing 
edge (ref. 16). A procedure similar to that used in reference 17 was employed. In each 
case an estimate was made of the point beyond which the inviscid velocity became vinrea- 
sonable because of the trailing-edge stagnation point. _ Then a least -square -curve fit of 
the velocities was made between that point (97-percent-chord station) and the 80- or 90- 
percent-chord station. This least -square -curve fit was then used to recompute the veloc- 
ities between the 80- or 90-percent station and the trailing edge. The outer surface was 
computed first so that the Nth value of the recomputed outer surface velocity (value closest 
to trailing edge) could also be included in the inner surface curve fit. A linear-curve fit 
was used on the outer surface and a parabolic fit on the inner surface. This procedure 
produced smooth and reasonable velocities for the boundary -layer calculation. It was 
determined that uncertainties in this method gave 2-percent or less error in /i. The fig- 
ures to be presented indicate the range in which the curve fit was made. 


It should be clearly understood, however, that this trailing-edge procedure is simply 
an approximate scheme in place of calculating the true effect of the correct wake, (For a 
nonzero trailing-edge thickness the total source strength is nonzero, and therefore some 
wake is simulated.) It may be possible to simulate the correct wake effect by use of a 
wake annular thickness distribution which has just enough geometric camber to cancel the 
AC due to thickness, thereby making the inner and outer pressures equal. 

Jr 


In order to obtain convergence of the boundary -layer displacement thickness, it was 
necessary to modify the wing coordinates, on the mth iteration, with the following effective 
displacement thickness 


6 = ^ ^ 


+ (m + 1)6 


m 


1 + 2 + 


+ (m + 1) 


(44) 


6. is based on velocities computed by using the actual wing geometry. 

6„/2, rather than 5^ as 


where 

first iteration the wing coordinates were modified by 
in equation (44). When separation occurred during the iteration process (shape 
H 2 2.4), the shape factor was limited to 2,55 to permit convergence. 


On the 

indicated 

factor 


Computer Program 

The complete computer program is described in appendix A. Input and output quan- 
tities and the program listing are included. One output quantity is the critical Mach num- 
ber, which is determined by an approximate method of von Karman (ref, 18), The printed 
value of the critical Mach number is valid only when M^ is low enough for the flow to be 
assumed incompressible. 
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The number of chordwise divisions N was 150. For the nacelles of this study the 
magnitude of the leading- edge suction peak had reached a stationary value by this value of 
N. For the mass-flow computation the radial increments were ^ijmer 
velocity generally had a very small variation across the flow area (substantiating the usual 
one-dimensional-flow assumption), these radial increments were more than sufficient. 

Calculations on a CDC-6600 computer required 2 minutes for an inviscid calculation 
and 4 minutes for a viscous calculation. Each additional iteration on the viscous calcula- 
tion required about 2 minutes. This time could be reduced by use of a cosine distribution 
of singularities which would increase the singularity density in regions of strong pressure 
gradients and reduce the density in the center of the body. Such a distribution of vortexes 
will produce the same downwash at points midway between the vortexes as a continuous 
distribution of vorticity (app. B of ref. 19). 

COMPARISON OF THEORY AND EXPERIMENT 

Calculations were compared with experimental data obtained in two wind tunnels 
(refs. 5 and 20) to provide a check on the present theory. The comparison with the tests 
of reference 5 is discussed first. These tests covered a range of subsonic Mach numbers 
at primarily zero angle of attack. Figure 3 (a) shows a sketch of the experimental setup 
and its relation to the theoretical model used to represent it. The mass flow in the experi- 
ment was controlled by the insertion of stationary sets of fan blades. The blockage of the 
sting, struts, and fan blades was represented in the theory by the blockage of a center 
body. Ih the experiment, the highest mass flow was obtained without a set of blades 
installed. In this condition the experimental setup was more faithfully modeled by the 
theory than with sets of blades installed, and the best agreement between measured and 
calculated internal pressures would be expected. This expected result was achieved as 
will be discussed. 

The pressure distribution and mass flow are not sensitive to variations in the diame- 
ter D of the vortex and source rings (see fig. 1). The value of D used in this study 
was the average of the minimum inside diameter and the maximum outside diameter of the 
wing. 

Nacelles 2 and 3 of reference 5 were selected for the present study (fig. 3(b)). These 
nacelles were designed to have prescribed external pressure distributions and internal 
contraction ratios. Comparisons between theory and experiment were made for three 
Mach numbers and several mass-flow ratios. The exact sting leading edge and nose shape 
were not given. Thus the nose of the center body in the theoretical model was placed at 
approximately the position of the sting leading edge and the shape was arbitrarily selected 
to be that of a NASA supercritical body of revolution (ref. 21). At the maximum diameter 
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Nacelle 




(b) Profile shapes for nacelles 2 and 3. 

Figure 3.- Geometry of experimental setup and theoretical model, 

station of the center body, a constant-diameter cylinder was attached and extended down- 
stream. As seen in figure 3(a), the constant diameter portion of the center body extended 
to infinity for the high-mass -flow cases of figures 4 to 14, However, for the low-mass- 
flow cases of figures 15 and 16 the center body was expanded at the approximate location 
of the fan blades. The necessity of this expansion for the low-mass -flow cases is discussed 
subsequently. 
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The quoted accuracy of the measured mass flows in reference 5 is 1.5 percent based 
on boundary layer and instrument errors. However, Young also states that the presence 
of the fan blades tended to separate the flow. Hence, for the lower mass flows tested, 
larger errors in the mass -flow ratio ju may be present. Errors in p would be reflected 
in a lack of agreement of measured and calculated pressures, especially on the inner sur- 
face and at the leading-edge suction peak. 

In the experiments the transition bands were applied from 6.25 to 10.4 percent of the 
chord. Therefore in the theory the boundary layer wa§ approximated by a turbulent bound- 
ary layer over the entire chord length, as discussed previously. The Reynolds number 
based on chord length was fixed at the experimental value of 2 x 10^. 

In all the pressure-distribution curves shown herein, with the exception of fig- 
ure 2(a), the upper and lower curves correspond to the outer and inner surfaces, 
respectively. 

A comparison of the theory and experiment is presented in figures 4 to 12. Mach 
numbers of 0,3, 0.5, and 0.7 were examined. In selecting the center-body diameter, two 
approaches are available. First, if the desired mass flow is known, the procedure is sim- 
ply to vary the body size until the desired mass flow is achieved. Only a few trial solu- 
tions are necessary to establish a curve of mass flow versus maximum body radius from 
which the maximum body radius can be selected to yield the desired mass flow. The 
second approach is to match the internal blockage area between the experimental geometry 
and a body of revolution and calculate the mass flow. Figures 4 to 7 used the first 
approach, and figures 8 to 11 used the second. Each figure gives the number of iterations 
on the boundary layer and the chordwise range from which velocities were used to obtain 
the least -square -curve fits. The iteration was continued until the maximum change in 6* 
on either surface was less than 2 percent of the maximum 6* for the respective surfaces. 
Good agreement is obtained between measured and calculated pressures when the boundary 
layer is included. 

Note that figures 6, 7, 10, and 11 have slightly supercritical flow as shown by the 
value of the critical pressure coefficient C* . A weak shock seems to occur near the 
leading edge. However, aft of this region agreement is good. 

Insist into the effect of the boundary layer on the mass flow is provided by exami- 
ning figures 8 to 11. These figures indicate that for a fixed geometry, the presence of the 
boundary layer increases the mass flow, contrary to what might at first be expected. 

Recall the case of a two-dimensional airfoil where the boundary layer causes an effective 
decrease in angle of attack; that is, the boundary layer reduces the lift and therefore the 
area under the curve of Cp versus 4 (ref. 17). As seen from figures 8 to 11, the 
boundary layer has the same effect on an annular wing, only in this case the radial force 
is being lowered. As this happens, the leading-edge stagnation point moves more toward 
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Figure 5.- Pressure distribution on nacelle 2; = 0.5; ^ 

Boundary layer: Six iterations; least square 90- to 97-percent c 




Figure 6.- Pressure distribution on nacelle 2; = 0.7; ju.^ = “ 

Boundary layer: Five iterations; least square 80- to 97-percent c. 



Figure 7.- Pressure distribution on nacelle 3; = 0.7; m+u = M = 

Boundary layer: Five iterations; least square 80- to 97-percent c 




Figure 8.- Pressure distribution on nacelle 3; M^o = 0.3; (^b/‘^)niax ~ percent; 

Pexp = 0.76. Boundary layer: Six iterations; least square 90- to 97-percent c. 
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x/c 


Figure 9.- Pressure distribution on nacelle 2; = 0.5; ( T,fc) = 12,61 percent; 

” \ /max 

f^exp “ Boundary layer: Six iterations; least square 90- to 97-percent c. 
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Figure 10.- Pressure distribution on nacelle 2; = 0.7; ( Ty.lc) = 12.61 percent; 

^ /max 

IXf^Tin = 0.74. Boundary layer: Five iterations, least square 80- to 97-percent c. 



Figure 11." Pressure distribution on nacelle 3; = 0.7; = 12.79 percent; 

Mexp = Boundary layer: Five iterations; least square 80- to 97-percent c. 
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O Nacelle 2 
A Nacelles 

Clear symbols: Without 
boundary layer 

Solid symbols: With 
boundary layer 


Figure 12,- Summary of calculated mass -flow ratios for cases where internal blockage 
was estimated, including effect of boundary layer. 


the outer surface. Thus far ahead of the body the capture streamline for the internal mass 
flow is at a larger radius and more flow is swallowed by the wing. A summary of the effect 
of the boundary layer on the mass flow is shown in the following sketch: 
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Both methods of selecting the center-body size show good results and neither seems 
superior to the other. The center-body size used in figures 8 to 11 was estimated from 
information given in reference 5. The accuracy with which the mass flow was calculated 
is summarized in figure 12 which shows that inclusion of the boundary layer in the calcu- 
lations improved this accuracy. 

Figures 13 and 14 compare the present theory with the methods of references 5 and 
10 for the flow condition of figure 4, Figure 13 is without the boundary layer and figure 14 
is with it. All methods are seen to give very close to the same results; therefore, the use 
of surface singularities gives essentially the same result in this case as the present 
method, which uses singularities placed on a mean diameter cylinder. 

When comparing different theoretical solutions of the flow over a nacelle, caution 
must be exercised to ensure that the mass flow from each solution is the same. Because 
the mass flow is a function of the circulation, different methods of approximating the Kutta 
condition can result in different mass flows on the same geometry. 

The experiments of figures 4 to 11 were run without fan blades present (highest 
mass -flow runs). When fan blades were present, it was found that the straight center- 
body model used previously was not sufficient to calculate accurately the internal nacelle 
surface pressures. Hence the center body was modified to approximate the experimental 
sting up to the fan-blade location, which was estimated to be at about the 7 5 -percent -chord 
location (the exact sting diameter and fan-blade location were not given in ref, 5). In this 
region the fan blade blockage was simulated by expanding the center body into an ellipsoid 
with a circular cross section (see fig. 3(a)), The ellipse in the x-r plane had a semi - 
minor axis of about 13 percent of the chord and a semimajor axis (r -direction), which was 
adjusted to give the experimental mass flow. Where the ellipsoid intersected the plane of 
the nozzle exit, a constant-diameter cylinder was extended to infinity. Figures 15 and 16 
show the results of using this model for two low-mass -flow cases. The external pressure 
distributions are well predicted in spite of the discrepancies on the internal surface. Fur- 
thermore a fairly reasonable estimate of the internal pressure distribution is obtained. In 
figure 16 some large experimental errors may occur in jix because of internal-flow sepa- 
ration as mentioned previously. This would explain the discrepancy at the leading-edge 
suction peak. 

Calculations were also compared with experimental data obtained from reference 20. 
The investigation of reference 20 was made on several NACA 1-series inlet contours tested 
over a range of subsonic and transonic Mach numbers, mass -flow ratios, and angles of 
attack. The thickness ratio (t/c)„_^ of the two nacelles of reference 20 selected for the 
present study was much smaller than the thickness ratio of nacelles 2 and 3 (2 percent com- 
pared with 11 percent). As shown in figure 3(a), the nacelles were sting supported with a 
throttle plug at the rear of the model for controlling the mass flow (the strut arrangement 
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Figure 13,- Pressure distribution on nacelle 3; 

11 . . = 0.76. All theories without boundary layer 





Figure 14.- Pressure distribution on nacelle 3; = 0,3 

fj.., = 0.76. Both theories with boundary layer. 



Figure 16.- Pressure distribution on nacelle 3; M = 0.30; fi,, = iJ. 

Lll 

Boundary layer: Five iterations; least square 90- to 97-percent c 



differed from that of fig. 3(a) in that fore and aft sets of struts were used). The plug was 
moved forward and backward to vary the mass flow. As previously, the blockage of the 
sting, struts, and plug was represented in the theory by the blockage of a center body. 
High-mass -flow ratios were selected for study to minimize the effect of the plug on the 
afterbody pressures. 

The models of reference 20 consisted of the same afterbody attached to various 
inlets. The inlets had NACA 1-series outer profiles and various inner profiles depending 
upon the contraction ratio. The two inlets selected for the present study were the NACA 
1-81-100 with a contraction ratio of 1.012 and the NACA 1-85-100 with a contraction ratio 
of 1.009. The numbers in the NACA designation indicate first the series, then 
in percent, and finally in percent. The quantity d is the highlight diameter 

less twice the inlet lip radius; X is the inlet length; and is the inlet (and nacelle, 

in this case) maximum diameter. 

The nose of the center body in the theoretical model was again placed at the position 
of the sting leading edge and had the same shape as used previously. (See fig. 3(a).) 

The tests of reference 20 were conducted without artificial boundary -layer transi- 
tion. A limited amount of data obtained with a transition band located 2.54 cm aft of the 
lip (the maximum nacelle diameter was 45.72 cm) showed no difference in pressure when 
compared with free transition data. Thus, it again seemed reasonable to assume a turbu- 
lent boundary layer over the entire chord length for the viscous calculations. The Reynolds 
number was a function of Mach number. For of 0.40 and 0.80 the Reynolds number 

was 7.437 x 10^/cm and 1.159 x 10^/cm, respectively. 

The comparison between theory and experiment is shown in figures 17 and 18 
for Moo of 0-80 and 0.40. The theoretical center-body size was selected to give the 
mass -flow ratio measured in the experiment. The data on the inner surface and the data 
labeled "other" were obtained from the author of reference 20. The data points labeled 
"other" are an average of measurements made at that x/c location on the sting, mass- 
flow rake, and inner surface (these measurements were not separately available). Good 
agreement was obtained on the outer surface for both cases when the boundary layer was 
included. The very large adverse pressure gradient on the inner surface of the inlet in 
figure 17 occurs at the higher Mach numbers. The boundary -layer calculation indicated 
flow separation in this region. Thus, the boundary layer on the outer surface did not con- 
verge to any better than a maximum change in 6* of 12 percent of the maximum 6*. 

This lack of convergence was due to erroneous changes in the thickness and camber caused 
primarily by incorrect values of 6* on the inner surface. Hence, flow separation and a 
possible influence of the mass-flow plug on the measured pressures at the trailing edge 
could explain the discrepancies between the theory and experiment in figure 17. 

In figure 18 the boundary layer converged on both surfaces to a maximum change in 
6* of less than 0.2 percent of the maximum 6* for the respective surfaces. The 
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Figure 17.- Pressure distribution on NACA nacelle 1-81-100; contraction ratio, 1.012; 
= 0.80; = ^exp “ Boundary layer: Three iterations; least square 

90- to 97-percent c. 
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Experiment 



Figure 18.- Pressure distribution on NACA nacelle 1-85-100; contraction ratio, 1.009; 
Moo ~ 0.40; = Mgjjp = 0.815. Boundary layer; Six iterations; least square 

90- to 97-percent c. 
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discrepancy at the leading-edge suction point could be due to inaccuracies in the leading- 
edge shape. The pressure in this region is very sensitive to slight errors in the manu- 
facture. The inner surface pressures are approximately predicted by the theory. The 
calculated velocities in this region are somewhat higher than the actual velocities because 
the center -body radius used in the theory was larger than the sting radius. 

Finally, figure 19 shows the accuracy of calculating the critical Mach number on 
nacelle 2 by using von Karman's method (ref. 18). The free-stream Mach number must 
be low enough to assume incompressible flow. The critical conditions are based on the 
external pressure distribution. The curve labeled "Minimum Cp" is from the leading- 
edge suction peak. The intersection of this curve and the curve for Cp* gives a critical 
Mach number for nacelle 2 of 0.632. The method of von Karm^ yields a critical Mach 
number of 0.658. 

In the case of a propulsive device the addition of energy to the internal flow will 
have an effect on the pressure distribution over the nacelle. Energy addition is treated in 
references 10 and 22. Shollenberger (ref. 22) has made a theoretical and experimental 
study of a two-dimensional propulsive lifting system. The external pressure distribution 
was found to have a relatively weak dependence on the intensity of energy addition to the 
internal flow; however, the inner surface pressure distribution was strongly affected. The 
theoretical formulation of the present study could be extended to include energy addition by 
use of a method similar to that of reference 22. 

CONCLUSIONS 

A method analogous to classical airfoil theory has been developed for computing the 
pressure distribution on an annular wing with a center body at zero angle of attack and 
subsonic speeds. The following conclusions are drawn: 

1. For the geometries investigated, this method was found to give essentially the 
same results as a method using a surface distribution of singularities and a stream-tube 
method. 

2. Comparison of the results with experimental data for a range of flow conditions 
and thickness ratios showed that when the boundary-layer displacement effects are included 
in the calculations, the pressure distribution and mass flow can be calculated with good 
accuracy. An accurate prediction is dependent on the existence of attached flow and a 
fairly accurate modeling of the internal-blockage geometry. 
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Figure 19.- Critical Mach number for nacelle 2 with fr, 



3. For cases in which the internal blockage was not accurately modeled or flow sepa- 
ration may have been present, the inner pressures were not accurately calculated. How- 
ever, the outer pressures were still calculated to a fairly good degree of accuracy. If 
internal flow separation occurs, the experimentally measured mass flow will be in error, 
and this error must be accounted for when comparing experimental and theoretical results. 

4. When calculating the boundary -layer thickness, the inviscid velocities were extrap- 
olated to the trailing edge. A more accurate calculation would include the effect of the 
wake. 

5. The method requires modest computer running times and storage and can be 
generalized to include such effects as lift and internal energy addition. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Hampton, Va., May 9, 1974. 
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APPENDIX A 


COMPUTER PROGRAM 
Input Variables 

The input data format is the Namelist method 'used by the Control Data series 6600 
computer system at the Langley Research Center. The input variables are listed in this 
appendix and the input format can easily be modified to any other format. The body radius 
and wing inner and outer coordinates are each listed in sequential order to correspond with 
the x-values listed. The x-values are input from the leading edge to the trailing edge. The 
first quantities in the arrays X, RB, RI, and RO are X(l), RB(1), RI(1), and RO(l), 
respectively. The positive direction for a radius is shown in figure 1; however x is 
measured from the leading edge, positive downstream. Any length units can be used so 
long as the same units are used for all variables. Meters are indicated, however, for 
dimensional quantities. 

The program uses 1.4 for y, the ratio of specific heats; however, a gas which has 
a value of y different from 1.4 can be treated by changing y in the pressure-coefficient 
equations and the mass -flow equation. Incompressible cases can be run by simply inputting 
zero for the free-stream Mach number. Special cases, such as camber with zero thickness 
or thickness with zero camber, can be treated by simply inputting the correct outer and 
inner surface coordinates. Some limitations are discussed. (See CAM.) For no center 
body, input all values of RB equal to 0 or set SCALE equal to 0 (both RB and SCALE 
must be input). 

AIAM 

C wing chord, m 

CAM use 1. when there is camber; use 0. when desire all Fq = 0, which is an 

inviscid thickness -only case with no center body. Note that zero geometric 
camber ^all o!o = 0^ is just a special case of equation (6). However, because 
of the method of defining the vortex strengths, would become inde- 

terminant if = 0. Thus the form of the equation for v^ will not 
permit use of equation (6) (as given by eqs. (8) or (12)) to accovint for the 
center-body interference for zero camber. Therefore a geometry which 
has a center body but no geometric camber ^ all Oq = o) cannot be treated 
by this computer program. Likewise, if portions of a cambered wing have 
a mean camber line with zero slope, a very small slope must be used in 
these regions. These limitations could be removed by replacing o/q7q by 
in the vortex equations ^see definition of y^^. 
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APPENDIX A - Continued 


DMAX 

DMIN 

FINI 

ITMAX 

M 

N 

NN 

NRHO 

RB 

RI 

RO 

SCALE 

UNRREF 

X 

XBN 


maximum wing outside diameter, m 
minimum wing inside diameter, m 

use 1. for the last data set and 0. for preceding data sets when more than one 
set in input; one data set is shown in sample input 

number of times boundary layer is computed; potential flow is computed 
ITMAX + 1 times; use 0. for inviscid solution 

free-stream Mach number 

number of wing panels in chordwise direction (same N used in text, see 
pages 10 and 17); maximum value is 170 

number of x-stations where quantities in arrays X, RI, RO, and RB are 
input; maximum value is 170 

number of annuli internal flow area is divided into to compute mass flow; max- 
imum value is 170 

array of center-body radii; for X £ XBN, use .0 for values of RB, m 
array of inner wing radii, m 
array of outer wing radii, m 

scale factor to expand or shrink elements of RB so as to vary mass flow; 
SCALE is ratio of desired body radius to input body radius 

Reynolds number per meter 

array of x-stations where arrays RI, RO, and RB are input, m 
x-station of center-body nose, m 
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APPENDIX A - Continued 


XAFT^ Between the x/ c values of XFOR and XAFT, a linear least-square-curve 
XFORJ fit is applied to outer surface velocities. Then the outer surface velocities 

are recomputed from XFOR to the trailing edge with this curve fit. Simi- 
larly, a parabolic least -square -curve fit is applied to the inner surface 
velocities between XFOR and XAFT and the Nth value of recomputed 
outer surface velocity. Then the inner surface velocities are recomputed 
from XFOR to the trailing edge with this curve fit. 

Sample Input 

The input for the viscous calculation of figure 4 is shown below. If more than one 
case is to be run, the complete data set must be repeated by using the proper values for 
FINI. The last card is identification information and uses format (6A10). 


3.6E0M M=150* MN = 5S» M=.T0, r = 4.4l01. 
XFOR=.QO* XAFT=.97. 

SCALE=1. 125. 

ITMAX=6, 

AIAM=1 ,2580. 

FINI=1., CAM=1., 

XBN=.8, 

USRREF = 453500. . 


x = 

0,0000. 

.0250, 

.0450, 

.0700 


,3528. 

.4410. 

.5292. 

.6174 


.8171. 

.8228, 

.8820. 

.970? 


1.4994, 

1 .5876, 

1.6758. 

1 .7640 


2.2933, 

2.4256, 

2.5138, 

?.6020 


3.3517, 

3.5281 . 

3 . 8604. 

3. 7927 


4.4101 , 





2.2093, 

2.1758, 

2.1 467, 

?.1170 


2.0073, 

1 .9888. 

1.9760 . 

1 .9698 


1.9797, 

1 .9802. 

1 .9852. 

1 .9927 


2.0095, 

2.0117, 

2.0140, 

2.0162 


2.0290. 

2.0321 . 

2.0338. 

2.0360 


2.0542, 

2.0586. 

2.0617. 

2.0648 


2.0749, 





2.2093, 

2.2428, 

2.2717, 

2.3010 


2.4019, 

2.4174, 

2.4302. 

2.4412 


2.4626, 

2.4631 , 

2.4682, 

2.47S7 


2.4995, 

2.5000, 

2.4991 , 

2.4969 


2.4602, 

2.4456, 

2.4355, 

2.4240 


2.3020, 

2.2671 , 

2.2389, 

2.210^ 


2.0785, 





0.0000, 

0.0000, 

0.8000. 

0.0000 


0.0000, 

0.0000, 

0.8000, 

0.0000 


.0315, 

.0358. 

.0620 . 

.0840 


.1345, 

.1390, 

.1415, 

. 1450 


.1500, 

.1500. 

.1500. 

.1500 


.1500, 

.1500, 

. 1 5 0 0 , 

.1500 


. 1500, 




YOUNGS COWL : 

1 




nMIN=3.93qs, 0MAX=8.. NRHO=170. 


.0882, 

.1200, 

.1500, 

. 1 764 , 

. 2646 , 

.7056, 

.7938, 

.8000, 

.8057, 

.8114, 

1 .0584, 

1 . 1466 , 

1.2348, 

1.3230, 

1.4112, 

1 .8522, 

1.9404, 

2.0286, 

2.1168, 

2.2051 . 

2.7343, 

2.8666, 

2.9548, 

3.0871 , 

3.2194, 

3.9691 , 

4. 1014, 

4.1896, 

4.2778, 

4.3660. 

2.1019, 

2.0822, 

2.0694, 

2.0603, 

2.0307, 

1.9711, 

1.9773, 

1 .9782, 

1 .9787, 

1.9792, 

1 .9985, 

2.0016, 

2.0038, 

2 . 0 056 , 

2 . 0073 , 

2.0184, 

2.0206, 

2.0228, 

2 . 0246 , 

2.0268. 

?.0396, 

2.0427, 

2.0449, 

2.0480, 

2.051 1 , 

2.0683, 

2.0710, 

2.0723, 

2.073^^, 

2.0749, 

2.3153, 

2.3344, 

2.3459, 

2.3537, 

2.3811, 

2.4509, 

2 . 46Q2 , 

2.4611 , 

2.4616, 

2.4621 , 

2.4823, 

2.4876, 

2.4920, 

2.4956. 

2.4982, 

P.4934, 

2.4889, 

2.4832, 

2*4761 , 

2 , 4686 , 

2.4054, 

2.3855, 

2.3714, 

2.3493, 

2.3263, 

2.1730, 

2.1443, 

2. 1P57, 

2.1067, 

2.0877, 

0.0000, 

0.0000, 

0.0000, 

o.nooo. 

0.0000, 

o.nooo. 

0.0000, 

0 .0000 , 

•0193, 

. 0 26 3 • 

.0980, 

.1080, 

.1160, 

.1230, 

.1290, 

. 1460, 

.1480, 

.1495, 

.1498, 

.1500, 

.1500, 

.1500, 

.1500, 

.1500, 

. 1 500, 

' . 1500, 

. 1500, 

.1500, 

. ] 500, 

. 1500, 
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APPENDIX A - Continued 
Output Vafiables 

Most of the input quantities are printed and identified with the symbols of the pre- 
ceding section. Other quantities have the meanings listed in this section. In the output 
for the wing surface solution, all perturbation velocities, including center -body velocities 
are for the equivalent incompressible flow and are on the wing surface ^ = 1 j . In the 

mass -flow computation output the perturbation velocities and the RAD are also for the 
equivalent incompressible flow. The pressure coefficients, total velocities, terms of 
equation (42), MFC and MU are for the compressible flow. Other output sections are 
either labeled as to compressibility or the nature is obvious (for example, GQ, etc. are 
obviously for the incompressible flow). The coordinate system of figure 1 is used except 
in the case of the center-body slopes as noted in the following list. Output which follows 
the boundary layer output is for the wing shape modified with the final 6* distribution. 


ALPHO 

“o 

at 

ALPHOO 

“o 

at 

ALPHQBB 

left-; 

ALPT 

“t 

at 

ALPTH 

“t 

at 

CPO 

s 

on 

CPI 

s 

on 

CPBA 


on 


’O] 


^bi^ 

03 


in second listing 


= 0] 


= 0 ]’ 


ference on the center body) 

CPBI isolated center body C ; this C would be theoretically correct if the wing 
were absent 


DDEL 

DELI 

DELO 


* * 

)„ - 
m m- 


6„ - 6„ _j, m 


6 from equation (44) on wing inner surface^, m 

6 from equation (44) on wing outer surface, m 
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APPENDIX A - Continued 


DELST “ 

DRDX dr/dx of center body at with x positive downstream 

DRDXH dr/dx of center body at with x positive downstream 
GAM 

GQ rt/Voo at 

GQB gq at 

H shape factor, 6 /e 

I i 

IT number of times, boundary layer has been computed; IT - 1 equals iteration 

number 

J i 

LAM X on equivalent incompressible wing 

MCR critical Mach number 

MFC where 

MU M 

RAD radii of annuli in mass -flow computation, m 

RIN RI - DELI (input values of RI used), m 

ROUT RO + DELO (input values of RO used), m 

RHOO 

RHOO-BODY at 
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APPENDIX A - Continued 


TEMPIO 

THE 

UOURO 

UOURI 

VTO 

VTI 

VIVIN 


VI/VIN 

WXIJUOO 

WXIJUOI 

WXIJUO 

WRIJUO 

WXIQJUO 

WXIQJUI 

WXIQJU 


quantity under summation sign of equation (42) at a given value of p^. 

9, m 

V/Vqo on wing outer surface at X-stations, used to compute boundary layer 

V/Vqj, on wing inner surface at X-stations, used to compute boundary layer; 
for x/c < 0.02, UOURI is set equal to UOURO(l) to obtain more reason- 
able values 

(v/v„)2 on outer surface of wing 

(v/v„)2 on inner surface of wing first time printed; Q V^o + u )Ky at 
p^. second time printed (square of velocity term under summation sign in 
eq. (42)) 

summation of equation (42) evaluated from j = 1 to the value of J listed for 
any given value of VIVIN. The summation progresses from wing inner sur- 
face toward center line, or center body if one is present, at fixed axial loca- 
tion. This sum is evaluated at = 0 (midchord) presently but this station 
can be changed by changing card numbers T170, T320, T650, and T870 in 
program MASFLOW. The last value of VIVIN equals Vl/VIN. 

V^o/^oo equation (42) 

^yl^oo on wing outer surface 

^yj^oo ^oj wing inner surface ^p^^ = ij 

summation in equation (13); for points where p^j ¥= 1 (required for computing 
mass flow), this quantity is ^.t ^qjjPqj 

v,,/V„ at 

at on wing outer surface ^p^. = 1^ 

at on wing inner surface ^p^. = 1^ 

same as WXIJUO only with respect to u. /V 

^9 yj 
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APPENDIX A - Continued 


WRIQJU at 


WXQJU 


’Voo 

at 

^oj’^ 

) . 
O] 


WRQJUO 


/Voo 

at 

^oj 

on wing outer surface 

("oi 

WRQJUI 

\i 

/V„ 

at 

^oj 

on wing inner surface 

{^oj 


WRQJU summation in equation (35); for points where p^j + 1 (required for mass 
flow), this quantity is at 

WXBQJU u^/V^ at 

WRBQJO Vb/V„ at 

XI {j 

XIO i„j 

XOJ/C ^Oj/® 


Program Structure and Listing 

This program is written in FORTRAN IV, Version 2.3 for the Control Data series 
6600 computer system with the run compiler and execution routines operating under 
SCOPE 3.0. The program requires 65,000g words of storage. 

The program operates in the OVERLAY mode. The main overlay program is 
called ANNULAR and has six primary overlay programs - GEOMTRY, MATRIX, 
PRESURE, BNDLYR, CRITM, and MASFLOW. The overlay structure and subprogram 
arrangement are shown in the following sketch: 
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APPENDIX A - Continued 


The function LOCF(A) used by subroutine MTLUP is a library function which 
returns the address of argument A. No listing is provided for LOCF. 

Numerous comment cards have been inserted in the main overlay to assist the user 
in understanding the flow of the computations. The purpose of each overlay and subpro- 
gram is described by comment cards at the beginning of the overlay or subprogram. 
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APPENDIX A - Continued 


OVERLAY(RING.0»O) A 10 

PROGRAM ANNULAR<INPUT»0UTPUT»TAPE5sINPUT*TAPE6=0UTPUT,TAPEl0) A ?0 

A 30 
A AO 

PROGRAM ANNULAR EMPLOYS ANNULAR WING THEORY TO COMPUTE THE A 50 

SUBSONIC FLOW OVER NACELLES INCLUDING THE COMPRESSIBILITY EFFECTS A 60 

And VISCOUS effects a to 

A BO 
A 90 

Any errors or problems encountered in using the program should be a 100 
Directed to dr. michael j. mann at nasa langley. area code 804-827 a 110 

“3711. A 1?0 

A 130 
A 140 

A CARD DECK AND DOCUMENTATION FOR THE PROGRAM ARE AVAILABLE FROM a 150 
COSMIC. UNIVERSITY OF GEORGIA. ATHENS. GEORGIA. 30601. A 160 

A 170 
A 180 

This program is written in coc fortran iv. version 2.3 . to run on a i9o 

CDC 6600 series COMPUTERS WITH THE SCOPE 3.0 OPERATING SYSTEM AND A 200 

LIBRARY TAPE. A 210 

A 220 
A ?30 

COMMON /CORD/ X ( 170 ) »R ( 170 ) .C.RAVG.NN.N.OMIN. T ( 1 70 ) .RB ( 170 ) .XBN A ?40 

COMMON /ONE/ ROUT (170 ) ,RIN ( 170 ) .M.OMAX .XI ( 1 70 ) .XIO ( 170 ) . ALPHO ( 1 70 ) A 250 

1.ALPHOO(170.1) .XBL(170) .ALPT(170) .ALPTH(170) »DRDX(170) .DROXhdTO) . A 260 

2LAM.IT.ITMAX A ?70 

COMMON /TWO/ IP.CAM.6AM(170).GQ(170).ALPOO(170,1) A 280 

COMMON /THREE/ VOUT(170) .VIN(170) A ?90 

COMMON /FOUR/ XOVC ( 170 ) .RO ( 170 ) .RI ( 170 ) .DELO < 1 70 ) .DEL I ( 1 70 ) , WLS ( 40 A TOO 

1) .DOLOOdTO) .DOLOI (170) .UNRREF.SAM.XFOR.XAFT A TlO 

COMMON /FIVE/ CPO(170).MCR A 320 

COMMON /SIX/ NRH0.RAD( 170) .VIVIN.RXX.RBXI A 330 

Dimension ident(6) a ^4o 

REAL LAM.KSQ.M.MFC.MCR. JOHN.MU A 350 

A 360 

Input a 370 

A 3H0 

Namelist /geom/ n.nn.m.c.omin.dmax.x.ro.ri.fini.cam.nrho.kb.xbn.sc a 3->o 

IAlE.UNRREF.ITMAX.AIAM.XFOR.XAFT a 400 

0 read (5.GEOH) A 410 

READ (5.100) IDENT A 420 

write (6.80) A 430 

Write ( 6.100 ) ioent a 440 

write (6.90) C.OMIN.OMAX.NRHO.UNRREF.ITMAX.N.M.AIAM.XFOR.XaFT. (X d A 450 

l).RO(I).RI(I).RB(I).I=l.NN) A 460 

Do 20 1=1.40 A 470 

WLS(1)=1. a 480 

0 Continue a 490 

Do 30 1 = 1 . NN A soo 

X0VC(I)=X(I)/C A SlO 

ROUT(I)=RO(I) A B?0 

RIN(I)=RI(I) A S30 

OELO(I)=0. A 540 

DELI(I)=0. a 550 

DOL0I(I)=0. A 560 

DOLOO(I)=0. A 570 

0 CONTINUE A ‘^BO 

A 590 

Scale center body a aoo 
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APPENDIX A - Continued 

C A 610 

Do 40 1=1 »NN A 620 

R8(I)=RB(I)*SCALE a 630 

40 CONTINUE A 640 

write (6*120) SCALE*(X(I)*RB(I)*I=1*NN) A 6S0 

C A 660 

C COMPUTE equivalent INCOMPRESSIBLE CENTER BODY A 670 

C A 680 

Do 50 1=1 »NN A 690 

RR(I)=RB(I)*S0RT<1.-M**2) a 700 

50 CONTINUE A 710 

IP=-1 A 720 

SaM=0. A 730 

IT=1 A 740 

C A 750 

C COMPUTE NACELLE MEAN CAMBER AND THICKNESS AND EQUIVALENT A 760 

c Incompressible shape - then compute nacelle camber and thickness a ?70 

c Slopes and center body slopes a 78 o 

C A 790 

60 CONTINUE A ^00 

Call overlay (4Hring*i,o*6hrecald a «io 

C A «20 

C COMPUTE CAMBER AND THICKNESS VORTEX STRENGTHS - 0AM (I) AND GO ( I ) A 830 

C A 840 

Call overlay (4HRING,2*0*6HRECaLD a 850 

C A 860 

C COMPUTE NACELLE AND CENTER BODY SURFACE PRESSURES A 870 

C A 880 

Call overlay (4hring*3,o*6hrecald a 89o 

C A ROO 

c tfst for inviscid only soln or end of boundary layer iteration a RlO 

C A <520 

Ir (IT.GT.ITMAX) go TO 70 A R30 

C A R40 

C EXTRAPOLATE FOR TRAILING EDGE VELOCITIES AND COMPUTE BOUNDARY A 950 

c Layer, modify nacelle coordinates for displacement thickness. a 960 

C A 970 

Call overlay (4Hring,4,o*6hrecall) a 980 

GO TO 60 A 990 

70 CONTINUE AlOOO 

C Aloio 

C CRITICAL MACH NUMBER A1020 

C A1030 

Call overlay (4HRIN6,5,0*6HRECaLL) A1040 

C A1050 

C COMPUTE VI/VIN BY NACELLE INTERNAL MASS FLOW INTEGRATION A1060 

C A1070 

CALL OVERLAY (4HRING,6*0*6HRECaLL') AIO 8 O 

C A1090 

C CALCULATE MFC AllOO 

C AlllO 

MFC=( (RXX**2-RBXI**2)/( (DMIN/2.)**2) )*VIVIN*( (1.*(.2*(M**2)*(1.-VI All20 
1VIN*»2) ) )**2.5) • A1130 

M|J=MFC/AIAM A1140 

Write (6*iio> mcr*vivin.mfc*mu auso 

C A1160 

c Check to see if there is another case auto 

C A1180 

IF (FINI.EQ.O.) go to 10 AU90 

STOP A1200 
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C A1210 

C Al?20 

80 FORMAT (IHl) Al?30 

90 Format (//16X»14HINPUT GE0METRY//4Xt2HC=*F9.3.4X*5H0MIN=,F9.3»4X,S A1?A0 

lHDMAX*»F9.3.AX.5HNRH0=»I9/4X.7HUMRREF=,EU.4»4X,6HlTMAX=,l3t9Xt2HN A1250 

2=* I3»4X»2HMs,F9,4»4X»5HAIAM**F9.4»4X»5HXF0R=.F9.4,AX»5HXAF T = »F9.4/ A1?60 

3/9X»1HX.8X»2HRO*8X.2HRI*8X,2HR8//<4F10,4) ) A1270 

100 Format <6aio> ai280 

no Format (1H0»AX»4HMCR=»F9.4»4X»7HVI/VINs»F9.4*4X»4HHFC=*F9.A*4X»3h' A1290 

IU=»F9.4//4X»56HMCR IS CORRECT ONLY IF M Is LOW ENOUGH hOR INCOMP a A 1300 

2SSUMP) A1310 

120 Format <1H1/2X,29HCENTER body expanded by SCALE*5X»6HSCALE=*F10.4/ A1320 

1///9X»1HX*8X,2HR8//(2F10.4) ) A1330 

End A1340- 
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APPENDIX A - Continued 


Function skk (z) b lo 

B 20 

Function skk computes the complete elliptical integral of the b 30 

First kind K(k> b ao 

B 50 

S, WAGNER» NASA AMES RESEARCH CENTER* MOFFETT FIELD CALIF. 9A035 B 60 

B 70 

REFERENCE C. HASTINGS* APPROXIMATIONS FOR DIGITAL COMPUTERS* P.172 B 80 

B 90 

Input z = k**2 b lOO 

B lio 

IF <Z) 10*20*30 B 120 

10 write (6*60) Z B 130 

STOP B lAO 

20 SKK=1.57079632679A896 8 ISO 

Return b i60 

30 E=1.0-Z B 170 

IF (E) 10*40*50 B 180 

40 write (6.70) Z 8 190 

STOP B 200 

50 SkK*( ( (E*.1451196212E-01*.3742563713E-01)*Et.3590092383E-0l J*E+.96 8 ?10 

166344259E-01)*E*1.38629436112-( ( ( (E*.441787012E-02*.3328355346E-01 B 220 

2>*E*.6880248576E-01)*E+.12498593597)*E*.5)*ALOG(E) B 230 

RETURN B 240 

C B 250 

C B 260 

60 FORMAT ( IHl * 10X.21HFUNCTION SKK K*»2 = *E16.9) B 270 

70 FORMAT ( IHl ♦ lOX *7HK»*2 = *E16.9* 13HK (K) INFINITE) B 280 

End 8 290- 


56 



APPENDIX A - Continued 


FUNCTION EKA <Z) C 10 

C C 20 

c Function eka computes the complete elliptical integral of the c 30 

c Second kino e(k> c ao 

C C 50 

C S, WAGNER* NASA AMES RESEARCH CENTER* MOFFETT FIELD CALIF. 9A03S C 60 

C C 70 

C reference C. HASTINGS* APPROXIMATIONS FOR DIGITAL COMPUTERS* P.175 C 80 

C C 90 

c Input z s k**2 c lOO 

C C 110 

IF IZ) . 10*20*30 C 120 

10 WRITE <6*60> Z C 130 

STOP C 140 

20 EKA»1. 570796326794896 C 150 

RETURN C 160 

30 Esl.O-Z C 170 

IF tE) 10*40*50 C 180 

40 EKAsl.O C 190 

RETURN C 200 

50 EKA=1.^( ( (.1736506451E-01*E*.4757383546E-01)*E*.6260601220e-01)*E* C 210 

1.44325141463)*E-( ( (.526449639E-02*E*.4069697526E-01)*E*. 9200180037 C 220 

2E-01)*E*.24998368310)*E*ALOG(E) C 230 

RETURN C 240 

C C 250 

C C 260 

60 FORMAT ( IHl * 10X*21HFUNCTION EKA K**2 * *E16.9) C 270 

End C 280- 
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SUBROUTINE MTLUP (X» Y»M»N»MAXtNTAB» I » VARI » VARO) 0 10 

C D 20 

c Multiple table lookcup subroutine, it computes y*f(x> from d 3o 

C A TABLE USING M»1 OR 2 FOR FIRST OR SECOND ORDER INTERPOLATION. 0 AO 

C D 50 

C MULTIPLE TABLE LOOKf-UP ON ONE INDEPENDENT VARIABLE TABLE D 60 

C D 70 

C USES AN EXTERNAL INTERVAL POINTER (I> TO START SEARCH D 80 

C I LESS THAN 0 WILLI CHECK MONOTONICITY 0 90 

C 0 100 

DIMENSION VARI(l). VAR0(MAX«1), Y(l). V(3), YY(2) 0 llO 

LOGICAL EX 0 120 

IF CM.EQ.O) GO TO 170 D 130 

IF IN.LE.l) GO TO 170 0 140 

Exs.F. 0 150 

IF (I.GE.O) GO TO 60 0 160 

IF (N.LT.2) GO TO 60 0 170 

C 0 180 

C MONOTONICITY CHECK D 190 

C 0 200 

IF <VAR1 (2)-VARI (1) ) 20.20,40 D 210 

C 0 220 

C ERROR IN MONOTONICITY 0 230 

C 0 240 

10 KsLOCFIVARKD) 0 250 

C 0 260 

C LOCF<X) IS A SYSTEM' ROUTINE THAT RETURNS THE ADDRESS OF X 0 270 

C 0 280 

Print 190 . j.k. <vari ( j) , j*i,N) 0 290 

STOP 0 300 

C 0 310 

C MONOTONIC DECREASING 0 320 

C D 330 

20 Do 30 Js2»N D 340 

IF (VARI (J>-VARI (J-I) ) 30.10.10 0 350 

30 CONTINUE 0 360 

GO TO 60 0 370 

C D 380 

C MONOTONIC INCREASING 0 390 

C 0 400 

40 DO 50 J=2.N 0 410 

If (VARI (J)-VARI (J-1) ) 10.10.50 0420 

50 CONTINUE D 430 

C D 440 

C INTERPOLATION D 450 

C 0 460 

60 IF (I.LE.O) 1=1 D 470 

IF (l.GE.N) I=N-l 0 480 

C 0 490 

C LOCATE I INTERVAL (X ( I ) .LE.X.LT.X ( I *1 ) > D 500 

C 0 510 

If ( (VARI(I)-X)*(VARI (I + D-X) ) 100.100.70 0 520 

C D 530 

C IN GIVES DIRECTION FOR SEARCH OF INTERVALS D 540 

C 0 550 

70 IN=SIGN(1.0,(VAR1(I»1)-VARI(I))*(X-VARI(1))) D 560 

C 0 570 

C IF X OUTSIDE ENDPOINTS, EXTRAPOLATE FROM END INTERVAL 0 580 

C D 590 

80 If ( (I»IN) .LE.O) GO TO 90 0 600 
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IF { (I+IN) .GE.N) GO TO 90 0 6l0 

IsI*IN 0 620 

If ( (VARI (I)-X)*(VARI (I + U-X) ) 100»100,80 0 630 

C D 6A0 

c Extrapolation d 6S0 

C D 660 

90 EX=.T. D <=.70 

100 If (M.EQ.2) GO TO 120 D 680 

C 0 690 

c First order d too 

C D 710 

Do 110 NT=1.NTAB D 720 

110 Y(NT)*(VARO<I*NT)*(VARl(I*l)-X)-VAROn+l»NT)*<VARI(I)-X))/(VARI(I+ 0 730 

11)-VARI (I.) > D 740 

IF (EX) I*I+IN D 750 

Return d t60 

C 0 770 

C SECOND ORDER 0 780 

C 0 790 

120 IF (N.EQ.2) GO TO 10 0 800 

IF (l.EQ.(N-l)) GO TO 140 0 810 

IF (I.EQ.l) GO TO 130 0 820 

C D 830 

C PICK THIRD POINT 0 840 

C 0 850 

SK=VARI (I + D-VARI (I) D 860 

IF ((SK«(X-VARI(I-l)>).LT.(SK*(VARI(I+2)-X))) GO TO 140 D 870 

130 Lsl D 880 

60 TO 150 0 890 

140 L=I-1 D 900 

150 V(1)=VARI(L)-X D 9l0 

V(2)»VARI (L+D-X 0 920 

V(3)*VARI (L+2)-X D 930 

DO 160 NT=1.NTAB ' D 940 

YY(1)=(VARD(L»NT)*V(2)-VARD(L*1»NT)*V(1) )/(VARI <L*1)-VARI (L) ) 0 950- 

YY(2)*(VARD<L*1»NT)*V<3)-VARD(L-^2,NT)*V(2) )/(VARI (L*2) -VAR I (L^ 1 ) ) D 960 
160 Y{NT)=(YY(1)*V(3)-YY(2)*V(1) )/(VARI(L*2)-VARI(L) ) 0970 

IF (EX) I=I+IN 0 980 

return D 990 

C DlOOO 

c Zero order oioio 

C 01020 

170 Do 180 NTS*1»NTAB D1030 

180 Y(NT)=VARD(1*NT) 01040 

Return oioso 

C D1060 

C 01070 

190 Format (ihi*49h table below out of order for mtlup at position »i 0108 O 

15,/31H X TABLE IS STORED IN LOCATION *06.// (8G15.8) ) 01090 

End 01100“ 
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OVERLAY(RIN6»1,0) E 10 

Program geomtry e 20 

C E 30 

c This overlay computes nacelle mean camber radii and thickness E 40 

C FOR THE EQUIVALENT INCOMPRESSIBLE SHAPE. THEN THE SLOPES ARE E 50 

C COMPUTED ON THE EQUIVALENT NACELLE CAMBER AND THICKNESS SHAPES E 51 

C and on the equivalent center body. E 52 

C E 60 

COMMON /CORO/ X ( 170 ) .R ( 170) »C»RAVG»NN.N.DMIN. T ( 170 ) »RB ( 170 ) .XRN E 70 

Common /one/ rout(17o) »rin(170) .m.omax.xi ( 170> »xio(i70 ) .alpho(17o> e bo 

1»ALPHOOU70*1) .X8Ln70) *ALPT(170).ALPTH(170) »DRDX(170) »DRDxH(170) , E 90 

2LAM.IT.ITMAX E lOO 

RFAL LAM.M E llO 

C E 111 

C COMPUTE NACELLE MEAN CAMBER RADII AND THICKNESS AND E ll2 

C COMPUTE NACELLE EQUIVALENT INCOMPRESSIBLE BODY E ll3 

C E 114 

DO 10 1=1. NN E 120 

R(I)=(ROUT(I)*RIN(I) )/2. E 130 

T(I)=POUT(I)-RIN(I) E 140 

10 CONTINUE E 150 

DO 20 1=1. NN E 160 

R(I)=R(I)»SQRT(1.-M**2) E 170 

T(I)=T(I)*SQRT(1.-M**2) E 180 

20 CONTINUE E 190 

DmAX=DMAX*SQRT(1.-M**2) E 200 

DMIN=OMIN*SQRT (l.-M**2) E ?10 

RAVG=(DMIN*DMAX)/4. E ?20 

C E ?30 

c Compute nacelle camber and thickness slopes and center body slopes e ?40 

c E 250 

CALL SLOPE (XI.XIO.ALPHO.ALPHOO) E ?60 

DO 30 J=l.N E 270 

XBL(J)=(C/2.)-(XI0( J)*RAVG) E 280 

30 CONTINUE E 290 

CALL THICK (ALPT.ALPTH) E 300 

CALL BODY (DRDX.DRDXH) E 310 

LaM=(DMAX*DMIN)/(2.*C) E 320 

IF (IT.LE.ITMAX) 60 TO 40 E 330 

Write (6.so) lam.(alpho<i).alphoo(I.1).xi(I).xio(I).alpt(I).alpth( e 340 

II) .DRDXm .DRDXH(I) .1.1 = 1. N) E 350 

40 Continue e 36O 

RETURN E 370 

C E 380 

C E 390 

50 FORMAT ( IHl .25X.42HGE0METRY ON EQUIVALENT INCOMPRESSIBLE BoDY//7x. E 400 

14hLAM=.F9.4//5X.5HALPH0.4X.6HALPH00.8X.2HXI.7X.3HXI0.6X.4HaLPT.5x. E 4l0 

25hALPTH.6X.4H0ROX.5X.5H0RDXH.4X.6HI OR J//(8F10.5.I10) ) E 420 

End E 430“ 
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Subroutine slope (xi«xio>alpho,alphoo) f lo 

c F 20 

c This subroutine computes xim, xio(j>, alphO(I)* alphoo(J) f 30 

C F AO 

C The XI AND XIO ARE FOR NACELLE AND CENTER BODY . F 50 

c The slopes are nacelle camber slopes. f so 

C F 70 

Common /cord/ xu 70) .Riiro) tc.RAVG.NN.NtOMiN f bo 

DIMENSION ALPHO(170)« ALPHOOI 170. 1 > » XK170)* X10(170) F 90 

lpa-1 F 100 

XXXs-,25*IC/N) F no 

XX=(C/N)*.7S F 120 

DO 30 I=1»N F 130 

If n.EQ.ll 60 TO 10 F lAO 

Call MTLUP (XXX»RXIi,2«NN*170»l«IP.X*R) F 150 

Call MTLUP (XX*RXI*2«NN»170tltlP»X.R> F 160 

TEMP4=(RXI1-RXI)/<C/N) F 170 

GO TO 20 F 180 

10 PlPssO. F 190 

P0P=<C/N)*.5 F 200 

call MTLUP <PIP.RXI1.2*NN*170*1*IP.X*R> F 210 

Call mtlup (pop.rxi.2,nn»17o»i,ip,x,r) f 220 

TEMP4s(RXI1-RXI)/( (C/N)*,5) F230 

20 CONTINUE F 240 

AlPH0(I)=ATAN(TEMP4) F 250 

XXX=XXX*(C/N) F 260 

XX=XX*(C/N) F 270 

30 CONTINUE F 280 

Do 50 1 = 1 »N F 290 

IF (I.EQ.N) 60 TO 40 F 300 

ALRH00(I»1)=(ALPH0(I)»ALPH0(I«1) )*.5 F 310 

GO TO 50 F 320 

40 AlPH00(I«1>=ALPH0(1)*(.5«(ALPH0(I)-ALPH0(I-1)>) F 330 

50 CONTINUE F 340 

TEMP5=.5*C/RAV6 F 350 

TeMP6=(C/N)/RAVG F 360 

Do 70 I=1»N F 370 

IF (I.EO.I) 60 TO 60 F 380 

XI(I)=XI(I-1)-TEMP6 F 390 

XI0(I)=X10(I-1)-TEMP6 F 400 

60 TO 70 F 410 

60 Xl (I)=TEMP5-(.25*TEMP6) F 420 

XI0<I)=TEMP5-(.75*TEMP6) F 430 

70 CONTINUE F 440 

RETURN F 450 

END F 460- 
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Subroutine thick <alpt,alpth> g lo 

6 20 

This subroutine computes nacelle thickness slopes, G 30 

G 40 

CONNON /com/ Xn70)»Rn70)»CtRAVG»NN.N'»DMIN»T(170) G 50 

dimension ALPTU70). ALPTH(170) G 60 

IP*-1 G 70 

XXX=-.25*<C/N) G 80 

XX®(C/N)*.75 G 90 

00 30 I=1»N G 100 

IF (I.EQ.l) GO TO 10 G 110 

Call htlup (xxx,txii,2,nn«17o.i<ip*x*T) g 120 

call NTLUP <XX*TXI*2«NN»170*1*IP*X*T> G 130 

TEMP4anTXIl-TXI)/2.)/(C/N) G 140 

GO TO 20 G ISO 

10 X14s(C/N)*,5 G 160 

Call MTLUP (Xl4»TXlt2»NNtl70»l,IP,X»T) G 170 

TEMP4s((0,-Txn/2,)/( (C/N)*,5) G 180 

20 CONTINUE G 190 

ALPT<I)aATAN(TEMP4) G 200 

XXXaXXX*(C/N) G 210 

XX=XX*<C/N) G 220 

30 CONTINUE G 230 

Do 50 laltN G 240 

IF (I.EQ.N) GO TO 40 G 250 

ALPTH(I)a<ALPT(I)*ALPT(I*l) )*.5 G 260 

GO TO 50 G ?70 

40 ALPTHa)=ALPT(I)*<.5*(ALPJ(I)-AtPT(I-l))> G280 

50 CONTINUE G 290 

Return g 300 

End 6 310- 
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SUBROUTINE BODY (DRDX.DRDXH) H 10 

C H 20 

c This subroutine computes center body slopes. h 3o 

C H 40 

common /cord/ XO70) .RdTO) »C»RAVG»NN»N»DMIN*T(170) .RB(170) tXBN H 50 

dimension ORDX(170)» DROXH(170) H 60 

lP=-l H 70 

XXX=-.25*(C/N) H 80 

XX=(C/N)*.75 H 90 

Do 30 I=l*N H 100 

IF (I.EQ.l) GO TO 10 H 110 

Call mtlup (xxx.rbxii»2»nn,17o,i.ip»x»rb) h 120 

IF <XXX.LE.XBN) RBXI1=0. H l30 

CALL MTLUP (XX,RBXI ♦2,nN* 170 . 1 » IP.XtRB) H 140 

IF (XX.LE.XBN) RBXI=0. H ISO 

DRDX(I)=(RBXI-RBXI1)/<C/N) H 160 

Go TO 20 H 170 

10 X14=(C/N)*.5 H 180 

Call mtlup (X14»RBX1»2.NN«170*1»IP»X»RB) H 190 

IF (X14.LE.XBN) RBXI=0. H 200 

DRDX(I)=R8X1/( (C/N)*,5) H2l0 

20 XXXsXXX*(C/N) H 220 

XX=XX*<C/N) H 230 

30 CONTINUE H 240 

DO 50 lal.N H 250 

IF (I.EQ.N) GO TO 40 H 260 

0r0XH<I)=(DR0X(I)+0RDX(I*1))/2. H 270 

Go TO 50 H 280 

40 DrDXH<I)=DRDX(I)*(.5*<DR0X(1)-DRDX(I-1))) H 290 

50 CONTINUE H 300 

IF (ORDXH(l).NE.O.) GO TO 80 H 310 

DO 60 I=2*N H 320 

IF (ORDXH(I) ,NE.O.) GO TO 70 H 330 

60 CONTINUE H 340 

70 ORDXH(I)sO. H 350 

80 CONTINUE H 360 

RETURN H 370 

End H 380- 
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OvERLAY<RIN6,2,0) I 10 

Program matrix i 20 

c I 30 

C This overlay computes nacelle camber and thickness vortex I 40 

C strengths - GAM(I) AND GO U ) I 41 

C I 50 

COMMON /CORD/ X ( 170 ) tR ( 170 ) »C»RAVG»NN.N*DMINtT ( 170 ) ,RB (1 70 ) *X8N ' I 60 

COMMON /ONE/ ROUT ( 170 ) ,R IN ( 1 70 ) *M ,DMAX ,X I { 170 ) ,XI0 < 170 ) t ALpHO < 170 ) I 70 

l.ALPHOO(170*l) »XBL(170) ,ALPT(170) ,ALPTH(170) *DROX(170 )*ORDxH(170) . I 80 

2LAM,IT.ITMAX I 90 

COMMON /TWO/ IP*CAMfGAM(170) .GQ{170) .ALPOO(170.1) I lOO 

Real ksq*lam*m i no 

Dimension wa(170>* wy(170)» gqb(170)» anscito). alphqbb{17o»1) . rh i 120 

100(170) I 130 

c I 131 

c Compute second alphoo which is printed - do 20 i 132 

c I 133 

DO 20 J*1.N I 140 

WRBQJUsO. I 150 

RHOO(J)=1. I 160 

XX=(C/N)*.25 I 170 

DO 10 I=ltN I 180 

Call mtlup (xx,rbxi,2,nn*17o.i,ip*x,rb> i i90 

IF (XX.LE.XBN) RBXI=o. I 200 

DELXIsXI (I)-XIO(J) I 210 

TeMP11=RH00(J)/(2.*( (0ELXI**2*RH00(J)**2)**1.5) ) I 220 

WRBQJU=WRBQJU+TEMP11*(RBXI/RAVG)*DRDX(I)*(XI(1)-XI (2) ) I 230 

Xx*XX*(C/N) I 240 

10 CONTINUE I 250 

AlPOO(J,1)sALPHOO(J»1) I ?60 

AlPH00(J*1)s-2.*3.14159*(ALPH00(J.1)+WRBQJU) I 270 

20 CONTINUE I 280 

C I 290 

C SOLVE matrix EQUATION FOR THE GAM ( I ) MATRIX I 300 

C I 310 

IF (CAM.EQ.O.) GO TO 50 I 320 

REWIND 10 I 330 

Do 40 J*1*N I 340 

Do 30 I=l»N I 350 

DELX1=XI (I)-XIO(J) I 360 

KSQ=4./(0ELXI**2*4.) I 370 

EK=EKA(KSQ) I 380 

SK=SKK(KSQ) I 390 

30 WA(I)=(((-1.)*DELXI*SQRT(KSQ))/4.)*<(((2.-KSQ)/(1.-KSQ))*EKA(KSQ)) I 400 

1-(2.*SKK(KSQ) ) )*ALPH0(I) I 4IO 

write (10) (WA(IT) .IT=1*N) »ALPHOO(J»l) I 420 

40 CONTINUE I 430 

Call ssleso (n«i«gam) i 440 

REWIND 10 I 450 

Go T& 7o I 460' 

50 DO 60 I=1*N I 470 

GaM(I)= 0 . 1 48 O 

60 CONTINUE I 490 

70 CONTINUE I 50 O 

C I 510 

C COMPUTE ALPHQBB - DO 80 I 52 O 

C I 530 

Do 80 Jsl*N I 540 

ALPHQBB (Jtl)sO, I 550 

XX=C/N*.2S 1 560 
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00 80 I=1»N I 570 

DELXIsXI(I)-XIO(J) I 580 

KSQs<»./(DELXI**2*4,) I 590 

EK=EKA(KSQ) I 600 

SK=SKK(KSQ) I 610 

Call MTLUP <XX»TX*2.NN,170»ltIP*X»T) I 620 

AlPHQBB(J*1)sALPHQBB(J«1>«( (TX/C)«( ( (SQRTil.-KSQJ »/2.)*{“0 eLX1/ABS I 630 

1 <DELXI) )»U(KSQ>2.)*SK)-<2.*<KSQ*1.>*EK) ) )• (XI ( 1 ) -XI (2) ) ) I 640 

XX=XX*(C/N) I 650 

80 CONTINUE I 660 

C I 670 

C SOLVE MATRIX EQUATION FOR GQB(D MATRIX I 680 

C I 690 

Rewind lo i ?oo 

Do 100 J=1»N I 710 

Do 90 I=1»N I 720 

DeLXI*XI(I)-XIO(J) I 730 

KSQ*4./(DELXI**2^4.> I 740 

EK=EKA(KSQ» I 750 

SK=SKK(KSD) I 760 

90 WY(I)=(((KSQ-2.)*EK>*2.*(1.-KSQ)*SK)/SQRT(KSQ))*((XI(1)-XI(2))/DELX I 770 

1I)*(1./LAM) I 780 

write (10) (WY(IT) .IT*1,N> .ALPhQBB(J*1) I 790 

100 CONTINUE I 600 

Call ssleso (n»i.gqB) i bio 

c I 620 

c Compute gq i b 30 

C I 840 

Xx=C/N*.25 I r50 

DO no I = 1«N I 860 

Call mtlup (xx.Tx.2»NNn70»i»ip.x,T) i bto 

G0(I)*((TX/C)/LAM)»(((C/(2.*RAVG))**2)*GQB(I)) I 880 

GQ(I>a-GQ(I) I 890 

XXsXX+(C/N) I 900 

no CONTINUE I 910 

IF (IT.LE.ITMAX) GO TO 130 I 920 

PRINT 140 I 930 

Do 120 I=1»N I 940 

write (6»150) I»GAM(I) ,00(1) tGOBd) .ALPHQBB(Itl) tALPHOOd*!) I 950 

120 CONTINUE I 960 

i3o' (Continue i 970 

RETURN I 980 

C I 990 

C IlOOO 

140 Format (lHl.llXdHI»9X,3H6AMdlX,2HGQ*10X,3HGQB.6X.7HALPHQBB»7X.6H IlOlO 

IAlPHOO) 11020 

150 FORMAT (9X. I3»2X.E12.5,4(2X»E11 .4) ) 11030 

End 11040- 
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SUBROUTINE SSLESO (NT»NCFLG» ANS) J 10 

c J ao 

c This subroutine solves simultaneous linear equations by j 30 

C SUCCESSIVE ORTHOGONALIZATION ** J AO 

C J 50 

C IF MATRIX ORDER IS * N, AND NUMBER OF RIGHT-HAND SIDES = M, THEN J 60 

C The amount of core that must be available for each array is — J 70 

C R * N^M, RV = N-»M'-1» CV = N-lt V =J(J^l) WHERE J IS THE INTEGER J 80 

C PORTION OF (N*M)/2 J 90 

C J 100 

c Dimensions based oni the values n*m ,le. 283t max n is 28i j no 

C J 120 

dimension RV<170>* R(170>» CV(l70)t V<8100)» ANS(NT) J 130 

REWIND 10 J lAO 

NlsNT*NCFLG J ISO 

J,=NI-1 J 160 

C J 170 

c Read (io> j i80 

c J 190 

READ <10) (R(I) .I*1.NI) J 200 

Do 10 I*ltJ J 2 l 0 

10 V(D»-R(I*1)/R(1) J 220 

IN=1 J 230 

C J 240 

C READ (10) J 250 

C J 260 

20 Read (10) (R(d »i®i»nd j ?.io 

12=0 J 280 

Do 40 I=1»J J 290 

RV(I)=0. J 300 

Do 30 I1=1»IN J 310 

12=12*1 J 320 

30 RV(I)=RV(I) ♦R(II)*V(I2) J 330 

N2=I*IN j 340 

40 RV(I)=RV(I) +R(N2) J 350 

I2=IN*1 J 360 

NN=J*IN*1 j 370 

KK=J*12 j 380 

J=J-1 J 390 

Do 60 Isl*J J 400 

DO 50 IIsl.IN J 410 

NN=NN-1 j 420 

KK=KK-1 j 430 

50 V(KK)=V(NN) j 440 

60 KK=KK-1 j 450 

DO 70 IsltlN J 460 

70 R(I)=V(I) J 470 

K=0 J 480 

Do 90 IsltJ J 490 

C=-RV(I*l)/RV(l) J 500 

DO 80 II=1«IN J 510 

CV(II)=C*R(II) J 520 

NN*K*II j 530 

12=12*1 J 540 

80 V(NN)=CV(II)*V(I2) J 550 

K=NN*1 j 560 

12=12*1 J 570 

90 V(K)=C J 580 

IF (J.EQ.NCFLG) GO TO 100 J 590 

IN=IN*1 J 600 
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Go TO 20 

J 

610 

100 

K=1-NT 

J 

620 


KK=0 

J 

630 


DO 110 Jsl*NCFLG 

J 

6^0 


K=K+NT 

J 

650 


kk=kk*nt 

J 

660 

110 

CONTINUE 

J 

670 


DO 120 J=1.NT 

J 

680 

120 

ANS(J)*-V(J) 

J 

690 


RETURN 

J 

700 


End 

J 

710 
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APPENDIX A — Continued 


OVERLAY<RIMG,3,0) K 10 

program presume K 20 

K 30 

This overlay computes nacelle and center body surface pressures K AO 

K 50 

COMMON /CORD/ X(170) »R{170) »C»RAVG»NN»N»DMIN»T(170) »RB<170) .XBN K 60 

COMMON /ONE/ ROUT ( 170) ,R1N ( 170) .M,DMAX,XI ( 170) .XIO ( 170) » ALPHO ( 1 70) K 70 

1»ALPHOO(170»1) ,XBL(170) ,ALPT(170) ,ALPTH(170) tORDX(170) »DRDxH(l70) , K 80 

2LaM.IT»ITMAX K 90 

COMMON /TWO/ IP»CAM»GAM(170) »GQ(170) »ALPOO(170*1) K lOO 

COMMON /THREE/ VOUT < 170) »VIN ( 170) K llO 

COMMON /FIVE/ CPO(I70),MCR K 120 

DIMENSION RHOO(I70), CPI (170) K 130 

REAL KSQtLAM.M K 140 

IF (IT.LE.ITMAX) GO TO 10 K 150 

PRINT 190 K 160 

Print i3o k i70 

PRINT 140 K 180 

PRINT 170 K 190 

0 CONTINUE K 200 

XXX=(C/N)*.75 K 210 

DO 120 J=1.N K 220 

XOJ=XIO(J)*RAVG/C K 230 

RH00(J)=1. K 240 

XX=(C/N)*,25 K 250 

K 260 

COMPUTE VELOCITIES ON NACELLE SURFACE - DO 20 K 270 

K 280 

WX1JU0=0. K 290 

WRlJUOsO. K 300 

WXlQJU=0. K 310 

WRlQJUsO. K 320 

WXQJU=0, K 330 

WROJU=0* K 340 

WXBOJU=0. K 350 

WRBOJUsO. K 360 

DO 20 I»1»N K 370 

DELXI*X1(I)-XI0(J) K 380 

KSQ=4./(DELX1**2*4.) K 390 

EKsEKA(KSQ) K 400 

SKsSKK(KSQ) K 410 

Call mtlup (xx,rbxi»2.nn»170,i,ip,x,rb) k 420 

IF (XX.LE.XBN) RBXI=0. K 430 

TEMP1=(1./(2.«SQRT(0ELXI**2*4.) ) )*( ( ( (2.-KSQ) / ( 1 .-KSQ) ) *EKa (KSQ) ) - K 440 

1 <2.*SKK(KSQ) )-( (KSQ/{1,-KSQ) )*EKA(KSQ) ) ) K 450 

TEMP2*( ( (-1,)*0ELXI*SQRT(KSQ) )/4.)*( ( ( (2,-KSQ) / ( 1 .-KSQ) ) *EkA<KSQ) ) K 460 

1-(2.*SKK(KSQ))) K 470 

TEMP7s(4,*EK)/( (SQRT(4,40ELXI**2) )*OELXI) K 480 

TeMP8=2,*TEMP1 K 490 

TeMP11=Rh6()(J)/(2.*U0ELXI**2*RH00(J)**2)**1.5) ) K 500 

TEMP12»(DELXI)/(2,*( (0ELXI**2*RH00(J)**2)**1.5) ) K 510 

Wx1JUOsWX1JUO>TEMP1*ALPHO(I)*GAM(I)*(1,/ <2, *3.14159)) K 520 

WR1JU0»WR1JU0^TEMP2*ALPH0(I)*GAM(I)*(1./(2.*3.14159) ) K 530 

WXlQJUaWXlQJU*(TEMPa*G0(I)*(XI(l)-XI(2))/(2.*3,14159)) K 540 

WRlQJUsWRlQJU*(TEMP2»GQ(I)*(XI(l)-XI(2))/(2.*3.14159)) K 550 

WXQJU=WXQJU-(TEMP7*ALPT(I)*(XI(1)-Xl(2))/(2.*3.14159)) K 560 

WRQJU=WRQJU*(TEMP8*ALPT(I)*(XI(1)«XI(2))/(2,*3,14I59)) K 570 

WXB0JUaWXB0JU*TEMPl2*(RBXI/RAVG)»DR0X(I)*(Xl(l)-XI(2)) K 580 

WRBQJU=WRBQJU*TEMPll*(RBXl/RAVG)*DROX<I)*(XI(l)-XI(2)) K 590 

XX=XX*(C/N) K 600 
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20 CONTINUE K 610 

IF U.EQ.N) 60 TO 30 K 620 

TENP3atALPMO<J)*6AM:<J>4ALPHO<J*l)*GAM(J*l) )*,25*<1,/(XI(1)-XI<2» ) ) K 630 
TEMP31»,25*(6Q<J)»GQ(J*1)) K 6<»0 

60 TO 60 K 650 

30 TEHP3sALPHO(J)*GAM(J)*.SO*tl./(XI(l)-XI (2)) >*.3333333333 K 660 

TeW»31=.50*<6Q<J>*. 3333333333) K 670 

60 Continue k 6S0 

WX1JUO0sWX1JU0«TEMP3 K 690 

WX1JUOISMX1JUO-TEHP3 K 700 

NX10JUI=WX1QJU-TEMP31 K 7l0 

WX10JUQstfX10JU«TEMP31 K 720 

WRQJU1sWRQJU«ALPTH(J> K 730 

WrQJUO=WRQJU-ALPTH(J) K 760 

C K 750 

C CALCULATE! TOTAL VELOCITY SQUARED AND APPLY COMPRESSIBILITY FACTOR K 760 

C K 770 

VT0*(l,*tWXlJU00*WXlQJU0*WXQJU»WXBQJU)/(l.-M**2) )**2 K 780 

VTI»(1,*<WXIJU0I*WX1QJUI*WXQJU*WXBQJU)/(1.-M**2) )**2 K 790 

C K 800 

C apply RIEGELS FACTOR K 810 

C K 820 

VTQ»VTO/ ( 1 . ♦ (T AN < ( ALPOO t J • 1) ♦ ALPTH ( J > ) /SORT ( 1 . -M**2 ) ) ) **2 ) K 830 

VTI*VTI/(1.*<TAN< (ALP00(J»1)-AlPTH(J) )/SQRT(1.-M**2) ) )**2) K 860 

VOUT(J)sSQRT(VTO) K 850 

ViN(J)=SQRT(VT1) K 860 

C K 870 

C COMPUTE CPO AND CPI K 880 

C K 890 

IF (M.EQ.0.) 60 TO 50 K 900 

CP0{J)s<2./(1.6*(M»*2) ))*< (!.♦(. 2*(M**2)*(l,-VTO) ) ) ** ( 1 .6/.6> -1 . ) K 9l0 
CPI<J)*(2./<1,6*(M**2) ) )*( (1,*(.2*(M**2)*(1.-VTI) ) )**a.6/.6)-l.) K 920 
60 TO 60 K 930 

50 CP0<J)=<1.-VT0) K 960 

CPI<J)a(l,-VTI) K 950 

60 CONTINUE K 960 

C K 970 

c Compute center body surface velocities - do to k 980 

C K 990 

WXlJaO. KlQOO 

WRlJaO. Kloio 

WXlQJsO. K1020 

WRlQJaO. K1030 

WXOJ-0. K1060 

WRQJsO. K1050 

WXBQJ^O. K1060 

WRBQJsO. KIOTO 

XX=(C/N)*.25 K1080 

Call mtlup (xxx»rbxj»2,nn»17o»i*ip»x»rb> kio9o 

IF IXXX.LE.XBN) RBXJ=0. KllOO 

RH 00 (J)»RBXJ/RAV 6 Kllio 

IF IRHOO(J) .LE.O.) GO TO 80 K1120 

DO 70 lal.N K1130 

DELXIaXI (I)-XIO(J) K1160 

Call mtlup <xx,rbxi,2»nn.170»i,ip,x,R8> kuso 

IF IXX.LE.XBN) RBXI=0. K1160 

KSQ«<6.*RHOO(J) )/( (0ELXI**2)*< (l.+RHOO(J) )**2) ) K1170 

EKsEKA(KSO) KllSO 

SKaSKK(KSQ) K1190 

TEMP1=(1./(2.*SQRT( (DELXI**2)+( (1,+RHOO(J) )**2) )))*(<( (2.-kSO)/(1. K1200 
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1 >KSQ ) ) *EK!) “ < 2 . »SK ) -( < 1 , /RHOO < J ) ) * ( KSO/ < I .-KSfl ) ) *EK ) ) K 1 2 1 0 

TEI«11P>2» C C 1-1 • ) ♦OELX 1 ) / ( 2 . *RHOO ( J> *SQRT < C OELX I **2 ) ♦ < U . ♦RHOO (J))**2) K1220 

ID }»U ( <2*-KSQ)/(l.-KS0) )*EK>-(2.*SK» ) K1230 

-TEMP7*(4,«0ELXI*EK)/( (SQRT( (DELXI**2>*( (l.*RHOO(J) )**2) > )*< <DELXI* K1240 

l*2)>ni.-RH00(Jn*»2) ) ) K1250 

TEHP8.*(2*/(RH00tJ)»SQRT((D£LXI**2)^Ul**RH00(J»>**2))))*(<(l..-<(2. K1260 

l*RHOO(J>*(RHOO<J)*l.) )/< (DELXI**2J ♦((l*-RHOO(J) )**2) ) ) )*EK)-SK) Kl2?0 

tEMPll=RHOOU)/<2.t( (DELXI**2>RH00(J)**2)**1.5) ) K1280 

TEMP12»(DELXl)/(2.t<(DELXl**2*RH00(J)**2)**l.5) ) K1290 

WX1J»WX1J*TEMP1*ALPH0 < I ) *GAM < I )* ( 1 ./ (2. *3, 14159) ) K1300 

WRl J=WR1J*TEMP2*ALPH0 ( I ) *6AM ( I ) * ( 1 ./ (2.*3, 14159) ) K1310 

WX1QJ*8X1QJ» (TEMP1*6Q(1)* (XI (l)-XI (2) )/(2, *3.14159)) K1320 

WR1QJ=WR1QJ*(TEMP2*6Q(1)*(XI<1)-XI (2) )/(2.*3.14159) ) K1330 

MXQJ=«XQJ- (TEMP7*ALPT ( 1 ) * (XI ( 1 )-X1 (2) ) / (2. *3. 14159) ) K1340 

WRQ J*WRQ J* ( TEMP8*ALPT ( I ) * ( X I ( 1 ) -X I ( 2 ) ) / ( 2 . *3 . 1 41 59 ) ) K 1 350 

WXBQJ*WXBQ J*TEMP 12* ( RBXI/RAV6 ) *ORDX ( I ) * ( X I ( 1 ) -X 1 ( 2) ) K 1 360 

WRBQJsWRBQJ^TEMPl 1* (RBXI/RAV6) *OROX ( I ) * (XI ( 1 ) -XI (2) ) K1370 

XX=XX^(C/N) K1390 

70 CONTINUE K1390 

80 VBA-(1.*(WX1J.WX1QJ^WXQJ.MXB0J)/(1.-H**2) )*»2^( (WRIJ.WRIQJ.NROJ^WR KI 4 OO 

1BQJ)/SQRT(1.-M**2) )**2 Kl4l0 

VBI=(1.^WXBQJ/(1,-M**2) )**2^(WRBQJ/SQRT(1,-M**2) )**2 K1420 

C K1430 

C COMPUTE CENTER BODY PRESSURES K1440 

C K1450 

IF (H.EQ.O.) GO TO 90 K1460 

CPBA=(2./(1.4*(M**2) ))*((!,♦(, 2*(M**2)*(1. -VBA) ) ) ** ( 1 ,4/,4) -1 . ) K1470 

CpBIs(2./(1.4*(M**2)) )*( (1,+(.2*(M**2)*(1,-VBI) ) )** ( 1 .4/,4) - 1 , ) KI 48 O 

GO TO 100 K1490 

90 CPBAs(l.-VBA) K1500 

CP8I*(1.-VB1) K151.0 

loO Continue kiszo 

IF (IT. LE. UMAX) go TO 110 K1S30 

WRITE (6.150) J.XOJ.CPO(J) tCPI (J) .VTO.VTI.WXIJUOO.WXIJUOI.WXIJUO.W K1540 

IRIJUO K1S50 

WRITE (6*160) WXIQJUO.WXIQJUI.WXIOJU.WRIQJU.WXQJU.WRQJUO.WRQJUI.WR K1S60 

IQJU K1S70 

Write ( 6 . 18 O) cpba.cpbi.wxboju«wrbqju.rhoo(J) kisso 

110 Continue kis9o 

XXXaXXX*(C/N) KI 6 OO 

120 continue K 1610 

return K1620 

C K1630 

C K1640 

130 FORMAT ( IX. 1HJ.7X.5HX0J/C.8X.3HCPO.9X.3HCPI .8X.3HVT0. 10X.3HVTI .6X. K1650 

17HWX1JU00.6X.7HWX1JUOI.7X.6HWXUU0.7X.6HWR1JUO) K1660 

140 FORMAT (19X.7HWX1QJU0.6X»7HWX1QJU1.7X.6HWX1QJU.7X.6HWR1QJU,8X.5HWX K1&70 

1 Q JU * 7X * 6HWRQ JUO . 7X . 6HWRQ JUI . 8X . 5HWRQ JU ) K 1 680 

150 FORMAT ( IX. 13.3 (2X.F10 ,5) .6 (2X.E1 1 .4) ) K1690 

160 FORMAT (13X.8(2X.E11.4) ) K1700 

170 FORMAT ( 19X.4HCPBA.9X.4HCPBI .7X.6HWXBQJU.7X.6HWRBQJU.4X.9HRH00-80D K1710 

IV/) K1720 

180 FORMAT ( 10X.2 ( 2X.F11 .5) .3X.2 (2X.E1 1 .4) ,Fl0.4/) K1730 

190 FORMAT ( IHl .20HNACELLE SURFACE SOLN. 10X.48HALL PERTURB VEL ARE ON K1740 

lEOUIV INCOMP NACELLE SURF.10X.26HCENTER BODY RHOO AND PRESS//) K1750 

End K1760- 
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APPENDIX A - Continued 


OVERLAY(RING*4,0) L 10 

PROGRAM BNDLYR L 20 

L 30 

This overlay extrapolates the velocities in the trailing edge L AO 

Region of the nacelle and computes the boundary layer on the l ai 

nacelle, then the nacelle coordinates are modified for the L A2 

DISPLACEMENT THICKNESS, L 43 

L 50 

common /CORD/ X(I70) ,R(170) ,C*RAVG»NN,N,DMIN,T(170) ,RB<170) ,XBN L 60 

COMMON /ONE/ ROUT ( 170 ) ,RIN (170 > . M.DMAX ,XI ( 170 ) .XIO ( 1 70 ) ♦ ALPHO (170 ) L 70 

ltALPHOO(l70.1) ,XBL(170) ,ALPT(170) ,ALPTH(170) .OROX(170) tOROxH(170) , L 80 

ELAfJf ITf ITMAX L 90 

COMMON /TWO/ IP.CAM.6AM(170) .GQ(170) ♦ALPOO(170,1) L 100 

COMMON /THREE/ VOUT ( 170 ) .VIN ( 170 ) L 110 

COMMON /FOUR/ XOVC ( 1 70 ) .RO ( 170 ) ,RI ( 170 ) .OELO ( 1 70 ) »DELI ( 170 ) , WLS (40 L 120 
1),DOLDO(170),DOLOI(170),UNRREF,SAM.XFOR»XAFT L 130 

DIMENSION UOURO(170), UOURK170)t THE(170), H(170), DELST(l70>» XL L 140 

1S(40)« YLS(40tl)> RESID(40»1), SUM ( 1) • ALS(3*3), BLS(3*1)« CLSO(40 L ISO 

2»2)» CLSI(40»3)» DOEL(170) L 160 

REAL JOHN»LAM»M L 170 

L 180 

EXTRAPOLATE OUTER SURFACE VELOCITIES FROM XFOR PERCENT CHORD TO L 190 

Trailing edge by linear least square fit of velocities l 200 

between XFOR AND XAFT PERCENT CHORD L 210 

L 220 

L=0 L 230 

DO 10 J*i*N L 240 

JOHN=XBL(J)/C L 250 

If (JOHN.lt. XFOR) GO TO 10 L 260 

If ( JOHN. GT, XAFT) 60 TO 10 L 270 

L=L>1 L 280 

XLS(L)sXBL(J) L 290 

YLS(L.1)=V0UT(J) L 300 

CONTINUE L 310 

Call LSQPOL (XLS*YLS»WLS*RESIO,L,SUM<l,ALS«BLSt2«CLSO,40,3) L 320 

DO 20 J=itN L 330 

JOHN*XBL(J)/C L 34() 

IF (JOHN. LT. XFOR) GO TO 20 L 350 

V0UT(J)»BLS(1*1)*BlS(2,1)«XBL(J) L 360 

CONTINUE L 370 

L 380 

EXTRAPOLATE INNER SURFACE VELOCITIES FROM XFOR PERCENT CHORD TO L 390 

TRAILING EDGE BY PARABOLIC LEAST SQUARE FIT OF VELOCITIES L 400 

Between xfor and xaft percent chord and last outer surface l 410 

VELOCITY FROM ABOVE L 420 

L 430 

LsO L 44O 

DO 30 Jsl*N L 450 

JOHN»XBL(J)/C L 460 

IF (JOHN. LT. XFOR) GO TO 30 L 470 

IF (JOHN, GT. XAFT) GO TO 30 L 480 

LsL.*l L 490 

xls(L):«xbl(J) l 500 

YLS(U»1)=VIN(J) L 510 

30 CONTINUE L 520 

LaL>l L 530 

XLS(U*X8L(N) L 540 

YLS(L»1)»V0UT(N) L 550 

CALL LSQPOL (XLS» YLS, WLS.RESIO.L .SUM* 1 , ALS.BLS* 3»CLSI »40f 3) L 560 

DO 40 Jxl.N L 570 
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J0HN»X8L(J)/C 



L S80 


IF (JOHN.LT.XFOR) GO TO 40 



L 590 


VlN(J)aBLS(l,l)^(BUS(2,l)*XBL( J) )+(8LS(3.1)*<X8L(J>**2) ) 


L 600 

40 

CONTINUE 



L 610 

C 




L 620 

c 

INTERPOLATE FOR VEL! AT INPUT X 



L 630 

c 




L 6<*0 


LUGsNN-l 



L 650 


DO 50 I=2tLUG 



L 660 


BlAIRsX(I) 



L 670 


CALL MTLUP (BLAIR»GLOSS»2*N»170»1»IP»XBL»VOUT) 



L 680 


Call mtlup (blair»tim»2»n*17o»i»ip»xbl»vin> 



L 690 


UOURO(I)=GLOSS 



L 700 


U0URI(I)=TIM 



L 710 

50 

CONTINUE 



L 720 


UOURO(l)sVOUT(l) 



L 730 


UOURI(I)=VIN(l) 



L 740 


UOURO(NN)=VOUT(N) 



L 7S0 


UOURI (NN)sVIN(N) 



L 760 


DO 60 I=I»NN 



L 770 


JOHN»X(I)/C 



L 780 


IF (JOHN. LT. .02) UOURI ( I ) =UOURO ( 1 ) 



L 790 

60 

continue 



L 800 

c 




L 810 

c 

COMPUTE OUTER BOUNDARY LAYER AND MODIFY OUTER WING 

coordinates 

FOR 

L 820 

c 

Displacement thickness 



L 830 

c 




L 840 


Call blyr (x,ro»nn,uouro.unrref»the.h»delst) 



L 850 


Do 100 1=1. nn 



L 860 


DEL0(I)*( (DEL0(I)*SAM)*(IT*DELST(I) ) )/(SAM*IT) 



L 870 


IF (IT.EQ.l) 70.80 



L 880 

70 

ROUT ( I ) =RO ( I ) ♦ (DELO ( I ) /2. ) 



L 890 


Go TO 90 



L 900 

80 

ROUT(I)=RO(I>*DELO(I) 



L 910 

90 

DOEL ( I ) =DELST ( I ) -DOLDO ( 1 ) 



L 920 


DOLDO(I)=DELST(I) 



L 930 

100 

CONTINUE 



L 940 

c 




L 950 

c 

Find new dmax 



L 960 

c 




L 970 


Do 140 Isl.NN 



L 980 


IF (I.EQ.1) 110.120 



L 990 

110 

ElIZ=ROUT(I) 



LlOOO 


GO TO 140 



LlOlO 

120 

IF (ROUT(I) .GT.ELIZ) 130.140 



L1020 

130 

ElIZ=ROUT(I) 



L1030 

140 

CONTINUE 



LI 040 


DMAX=2.*ELIZ 



LI 050 


write (6,230) IT. DMAX, (XOVC ( I ) ,X ( I )-*UOURO ( I ) .DELST ( I ) ,DDEL ( I ) , 

DELO 

L1060 


1 (l).ROUT(I).THE(I).H(I).I = l.NN) 



L1070 

c 




L1080 

c 

COMPUTE INNER BOUNDARY LAYER AND MODIFY INNER WING 

coordinates 

FOR 

L1090 

c 

DISPLACEMENT THICKNESS 



LllOO 

c 




LlllO 


Call blyr (x.ri.nn.uouri.unrref.the.h.delst) 



L1120 


Do 180 Isl.NN 



L1130 


deli (1 ) = ( (DEL I ( I ) *SAM ) ♦ ( IT*DELST (!)))/( SAM* IT ) 



L1140 


IF (IT.EQ.l) 150.160 



L1150 

150 

RlN(I)=RI(I)-(DELl(I)/2.) 



L1160 


GO TO 170 



L1170 
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160 RlN(l)sRIfl)-OCLl(I) LllSO 

170 OOELmsOELSTm-OOLOim L1190 

OOLOimsOELST(I) L1200 

180 CONTINUE L1210 

C LI 220 

C FIND NEW ONIN LI 230 

C LI 240 

00 220 Isl*NN LI 250 

IF (I.EQ.ll 190*200 L1260 

190 ELlZ-RINfll L1270 

60 TO 220 LI 280 

200 IF CRINID.LT.ELIZ) 210*220 L1290 

210 ELIZsRINd) L1300 

220 Continue li310 

0mIN»2.*ELIZ L1320 

write (6*240) IT*DH'IN*(X0VC(I)*X(I)*U0URI(I)*0ELST(I)*0DEL(I)*0ELI L1330 
1(I)*RIN(I)*THE(I)*H(I)*Is1*NN) L1340 

SAMsSAN^IT L1350 

lTsIT*l L1360 

Return l1370 

C L1380 

C L1390 

230 Format (1H1*4X*44H0UTER wing shape modified by boundary layer. /4X* L1400 

ISOHVELOCITIES ARE ONES USED TO COMPUTE BOUNDARY LAYER//4X*3H1T=* 13 Ll4lO 

2*4X.5H0MAX=*F9.4//7X*3HX/C*9X*1HX*5X*5HU0UR0*5X*5HDELST*6X,4H00EL* L1420 

36x*4HDELO*6X*4HROUT*7X*3HTHE*9x*IHH//( 9F10.4) ) Ll430 

240 FORMAT (IHl *4X*44HINNER WING SHAPE MODIFIED BY BOUNDARY LAyER./4x. LI440 

ISOHVELOCITIES ARE ONES USEO TO COMPUTE BOUNDARY LAYER//4X*3HIT«* 13 LU50 

2*4X*5HDMIN=*F9,4//7X*3HX/C*9X*1HX*SX*5HU0URI*SX*5H0ELST*6X,4H0DEL* L1460 

36x*4HDELI*7X*3HRIN*7X*3HTHE*9X*1hH//<9F10.4) ) L1470 

End L1480 
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Subroutine blyr (S,r»nx»uour»unrref»the»h,delst) m lo 

C H 20' 

c This subroutine compu(ES the incompressible turbulent axisymmetric m 30 

C BOUNDARY LAYER USING A SLiIGHT MODIFICATION OF TRUCKENBRODTs METHOD M 40 

c Page sba - sbs of boundary layer theory by schlichting»4Th edition m so 

C M 60 

C TRANSITION IS ASSUMED AT THE LEADING EDGE. M 70 

C M 80 

REAL lARGf IARGL*IU0URYL M 90 

DIMENSION R(200)« S(200) * UOUR(200) * IUOUR(200)« UAOUR(200)< REY(2 M lOO 

100>f CF(200>* AR6(200>« IARG(200). THE(200)» RTHE(200)» B<200>t AR M IlO 

26l( 200)» IAR6L(200)» L(200>« TU<13)t TH(13)* H(200)« OELST(200)« S M 120 

3QRTXI(200) M 130 

Do 10 Isl»NX M 140 

IF (UOUR(I) .LT.. 00001) UOUR ( I ) sO .00001 M ISO 

10 CONTINUE M 160 

C M 170 

C U AVE OVER U REF M 180 

C M 1^0 

Call INT (S*UOURtNX*IUOUR) M 200 

UAOUR ( 1 ) sUOUR ( 1 ) M 210 

Do 20 Is 2 «NX M 220 

20 UAOUR(I)=IUOUR(1)/S(I) M 230 

C M 240 

C CF M 2S0 

C M 260 

REY(1)sO. M 270 

CF(l)»l. M 280 

DO 30 Is2«NX M 290 

REY(I)=UNRREF*UA0UR>(I)#S(I) M 300 

IF (REY(I) .LE.l.) REY(I)»1. M 310 

CF(I)».0745/(REY(I)**.218)*. 00072 M 320 

30 CONTINUE M 330 

C M 340 

c Theta m 350 

C M 360 

DO 40 1*1. NX M 370 

40 AR6(I)*UOUR(I)**(10./3.)*R(I)**(7./6.) M 38 O 

Call int (s.arg.nx.iarg) m 390 

ThE(1)*0. M 400 

DO so I*2.NX M 410 

IF (lARG(I) .LE.O.) IARG(I)*0. M 420 

50 THE(I)*(UAOUR(I)*S(lj )**(!. /7.>*CF(I)/2.*IARG(I)**<6./7.)/(U0UR(I) M 430 

1**3*R(1) ) M 440 

C M 450 

c Shape factor m 460 

C M 470 

C XI M 480 

C _M 490 

SQRTxi(i)*o. m‘ ido* 

DO 60 I*1.NX M 510 

S0RTXKI)*(IARG(I) )»*2 M 520 

60 CONTINUE M 530 

C M 540 

CL M 550 

C M 560 

DO 70 1=1. NX M 570 

RTHE(I)=UNRREF*U0UR'(1)*THE(1) M 580 

IF (RTHE(l) .LT.l.) RTHE(I)=1. M 590 

B(I)=.07*AL0610(RThE<I) )-.23 M 600 
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70 


BO 

90 

100 


no 


Ir (B(I),LT. 0 .) B(I)= 0 , 

ArGL(I)=<B(I)-ALOO(UOUR(I) ) )*2.*SQRTXI (I) 

Call int (sqrtxi»argl»nx»iar6d 

L(1)=0. 

Do 100 I=2«NX 

IF ISQRTXI (I) .EQ.O.) GO TO BO 
L(I)=ALOG(UOUR(I) ) ♦IARGL(I)/SQRTXI(I)**2 
GO TO 90 
L<n= 0 . 

IF (Lm.GT.l.) L(I)=1. 

If (HI) .LT.-.19) LI<I)=-.19 

H AND DELTA STAR 

Data TL/-.133.-.172,-.157»-.131»-.1»-,05.,0».05».1*.3».5,.7»1,/.TH 
1/2.A.2.2*2.0.1.8»1.6A5,1.A95»1.4,1.335»1.295*1,2*1.15»1.12,1.07/,N 
2TL/13/ 

IP=-1 

DO no 1 = 1. NX 
TeP1=L(I) 

Call mtlup <tepi,tep2,2,ntl»13,i,ip,tl.th) 

H(I)=TEP2 

DELST(I)=THE(I)*H(I» 

CONTINUE 

RETURN 

End 


M Gio 

M A20 
M 630 
M 6A0 
M 650 
M 660 
M 670 
M 680 
M 690 
M 700 
M 710 
M 720 
M 730 
M 7A0 
M 750 
M 760 
M 770 
M 780 
M 790 
M BOO 
M alo 

M 820 
M 830 
M 8A0 
M 850 
M 860- 
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SUBROUTINE INT (X*Y»N»ANS) 


N 

10 

c 



N 

20 

c 

This subroutine is used by subroutine 

blyr to integrate y ox from 

N 

30 

c 

X(l> TO X(I>. RESULT EQUALS ANSd). I 

RUNS FROM 1 TO N. 

N 

^0 

c 



H 

SO 


Dimension xuoo)» y(ioo)* ans(ioo) 


N 

60 


SUM=0,0 


N 

70 


AnS(1)=0. 


N 

80 


If (n.eq.Z) go to 90 


N 

90 


Isl 


N 

ioo 


GO TO 40 


N 

110 

10 

XAsX(I)-X(I-l) 


N 

120 


XB*X(I*1)-X(I) 


N 

130 


Ya=Y(I)-Y(I-1) 


N 

140 


YBsY(I»1)-Y(I) 


N 

150 


IF (XA.EQ.O..OR.XB.EQ.O.) 60 TO'20 


N 

160 


IF (YA/XA.GE.2.*YB/XB*AND.YA/XA.GE.0. 

) GO TO 20 

N 

170 


IF <YA/XA.LE.2.*YB/XB.AND.YA/XA.LE.0. 

) GO TO 20 

N 

180 


SUMFW=XB*(Y(I) ♦YB/3.*(XA*»2*YB*XB**2*YA)/(6.*XA*(XA4XB) » > 

N 

190 


GO TO 30 


N 

200 

20 

SUMFW=XB*( (Y(I>*Y(I*1) )/2.) 


N 

210 

30 

IF <I.EQ.N-I) GO TO' 80 


N 

220 

M) 

Isl^l 


N 

230 


Xa=X(I)-X(I-1) 


N 

240 


XB=X(I*1)-X(I) 


N 

250 


Ya*Y(I)-Y(I-1) 


N 

260 


YB=Y(I*1)-Y(I) 


N 

270 


IF (XA.EQ.O..OR.XB.EQ.O.> 60 TO 50 


N 

280 


If (YB/XB.GE.2.«YA/XA.AND.YB/XB.GE.0. 

> GO TO 50 

N 

290 


If (YB/XB.L£.2.*YA/XA.AND.YB/XB.LE.0. 

) GO TO 50 

N 

300 


SUMBKsXA* (Y ( I ) -YA/3.- (XB**2*YA4XA**2*YB» / (6.*XB* (XA+XB) ) ) 

N 

310 


Go TO 60 


N 

320 

50 

SUMBK=XA*< (Y(I)»Y(I-l) )/2.» 


N 

330 

60 

IF (I.EQ.2) GO TO 70 


N 

340 


SUM=SUM*(SUMFW+SUMBK)/Z. 


N 

350 


AnS(I)=SUM 


N 

360 


60 TO 10 


N 

370 

70 

SUM»SUMBK 


N 

380 


AnS(2)sSUM 


N 

390 


60 TO 10 


N 

400 

BO 

SuMsSUH^SUMFW 


N 

410 


ANS(N)aSUM 


N 

420 


60 TO 100 


N 

430 

90 

AnS( 2)*( (Y<1)*Y<2) )/2.)*(X(2)-X(l) ) 


N 

440 

100 

CONTINUE 


N 

450 


RETURN 


N 

460 


End 


N 

470 
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SUBROUTINE LSQPOL (X«Y*WtRESIOyN«SUHfLt AfBtM«C*NMAX«MMAX) 0 10 

0 20 

This subroutine is a least square polynomial pit. given a set of O 30 

values of an independent variable X WITH ASSOCIATED WEIGHTS W AND 0 40 

A SET OF CORRESPONDING VALUES OF Y. THE ROUTINE DETERMINES THE 0 50 

COEFFICIENTS OF THE POLYNOMIAL OF DEGREE M-l WHICH GIVES ThE BEST 0 60 

FIT IN THE LEAST SQUARES SENSE TO THE SET OF Y, 0 70 

0 80 

DIMENSION X(NMAX)* Y(NMAXfU* RESIO<NMAX«L) • A (MMAX*MMAX) « B(MMAX« 0 90 

1L)» C(NMAX»M)y SUM(D* W(NMAX) 0 lOO 

10 DO 20 I»lfN 0 110 

20 C(I»1)*1,0 0 120 

Do 30 Js 2 ,M 0 130 

DO 30 IsltN 0 140 

30 C(I»J)sC(I»J-l)*X<I> 0 150 

DO 40 IrltM 0 160 

Do 40 JalfM 0 170 

A(I»J)sO.O 0 180 

Do 40 K=1»N 0 190 

40 A(I,J)sA(I»J)+C(K.I)«C(K»J)*W(K) 0 200 

Do 50 Jsl»L 0 210 

Do 50 I*1»M 0 220 

B(I*J)»0.0 0 230 

DO 50 Ksl.N 0 240 

50 B(I»J)sB(I.J)*C(K»I)*Y<KtJ)»W(K) 0 250 

call MATINV (A,M,B»L»DETERM.RESlD.CfMMAX.ISCALE) 0 260 

DO 70 J=1»L 0 270 

SUM (J) =0.0 0 280 

KK=M 0 290 

Do 60 K=1*M 0 300 

C(K»1)»B(KK, J) 0 310 

60 KK=KK-1 0 320 

Do 70 l=l*N 0 330 

RESID(I,J)=P0LYE1(X(I),M»C)-Y<I.J) 0 340 

70 SUM(4)=SUM(J)*RESID<I*J)**2*W(I) 0 350 

RETURN 0 360 

End 0 370- 
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SUBROUTINE MATINV (A»N,B*Nf DETERM, IPIVOT* INDEXtNMAX* ISCALE) P 10 

C P 20 

C THIS SUBROUTINE SOLVES THE MATRIX EQUATION AX*B WHERE A IS A P 30 

c Square coefficient matrix and b is a matrix of constant vectors. p ao 

c P 50 

dimension IPIVOT(N), A(NMAX,N), B(NMAXfM), INDEX (NHAX,2) p 80 

EQUIVALENCE (IROWfjROW), < ICOLUM, JCOLUM) , ( AM AX, T, SWAP) P 90 

C p 100 

C INITIALIZATION p HO 

C P 120 

10 ISCALE^O P 130 

Rl«10. 0**100 p 140 

R2«1.0/R1 p ISO 

Determ*!, 0 p i6o 

DO 20 J*1,N P 170 

20 IPIVOT(J)»0 p 180 

Do 360 1*1, N p 190 

C P 200 

C SEARCH FOR PIVOT ELEMENT P 210 

C P 220 

AmAX*0.0 P 230 

Do 70 Jsl,N p 240 

IF (IPIVOT(J)-I) 30,70,30 p 250 

30 Do 60 K*1,N p 260 

If (IPlVOT(K)-l) 40,60,400 p 270 

40 Ir (ABS(AMAX)-ABS(A(J,K> ) ) 50,60,60 p 280 

50 lROW*J p 290 

ICOLUMaK p 300 

AmAX*A(J,K) P 310 

60 CONTINUE p 320 

70 CONTINUE p 33O 

If (AMAX) 90,80,90 p 34O 

80 OETERMaO.O p 350 

ISCALEaO p 360 

GO TO 400 p 370 

90 IPIV0T(1C0LUM)*IPIV0T(IC0LUM)*1 p 38O 

C p 390 

C interchange ROWS TO PUT PIVOT ELEMENT ON DIAGONAL P 400 

C P 410 

IF (IROW-ICOLUM) 100,140,100 p 420 

100 DETERMa-DETERM p 43O 

DO no L*1,N P 440 

SWAPaA(IROW,L) P 450 

A(IROW,L)aA(ICOLUM,L) P 46O 

no A(ICOLUM,L)aSWAP p 47O 

IF <M) 140,140,120 p 480 

120 Do 130 L=1,M p 490 

SWAP*B( IROW,L> P 500 

8{IR0W,L)aB{IC0LUM,L) p 5IO 

130 B ( ICOLUM, DaSWAP p 520 

140 lNDEX(l,l)aiROW p 530 

lNOEX(I,2)alcOLUM p 540 

PlVOTaA( ICOLUM, ICOLUM) P 55O 

IF <PIVOT) 150,80,150 p 560 

C p 570 

C SCALE THE DETERMINANT p 58O 

C p 590 

150 PiVOTIaPIVOT p ftOO 
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IF <ABS<DETERM)-R1) 180*160*160 

P 610 

160 

0ETERM=0ETERM/R1 

P 620 


ISCALE=ISCALE*1 

P 630 


IF (ABS(OETERM)-Rl) 210*170*170 

P 640 

170 

DETERM=DETERM/R1 

P 650 


ISCALE*ISCALE+1 

P 660 


GO TO 210 

P 670 

160 

If <ABS(DETERM)-R2) 190*190*210 

P 680 

190 

DETERM=OETERM*Rl 

P 690 


ISCALE»ISCALE-1 

P 700 


IF <ABS(DETERM)-R2) 200*200*210 

P 710 

200 

oeterm=oeterm*ri 

P 720 


ISCALE*ISCALE-1 

P 730 

210 

If (ABS(OIVOTI)-Rl) 240*220*220 

P 740 

220 

PlVOTI»PIVOTI/Rl 

P 750 


I SCALE* I SCALE* 1 

P 760 


IF (ABS(PIVOTI)-Rl) 270*230*230 

P 770 

230 

PIVOTI*PIVOTI/Rl 

P 780 


ISCALE=ISCALE*1 

P 790 


Go TO 270 

P 800 

240 

IF (ABS(PIVOTI)-R2) 250*250*270 

P fllO 

250 

PIV0TI»PIV0TI*R1 

P 820 


ISCALE=ISCALE-1 

P 830 


IF <ABS(PIV0TI)-R2) 260*260*270 

P 840 

260 

PIV0TI=PIV0TI*R1 

P 850 


ISCALE*ISCALE-1 

P 860 

270 

determ*oeterm*pivoti 

P 870 

C 


P 880 

c 

Divide pivot row by pivot element 

P 890 

c 


P 900 


A(ICOLUM*ICOLUM)*1.0 

P 910 


DO 280 L*1*N 

P 920 

2B0 

A(ICOLUM*U=A(ICOLUM*L)/PIVOT 

P 930 


IF (M) 310*310*290 

P 940 

290 

Do 300 L=1*M 

P 950 

300 

B(ICOLUM,L)=B(ICOLUM*L) /PIVOT 

P 960 

C 


P 970 

C 

REDUCE NON-PIVOT ROWS 

P 980 

C 


P 990" 

310 

Do 360 L1*1*N 

PlOOO 


IF <L1-IC0LUM> 320*360*320 

PlOlO 

320 

T=A(L1*IC0LUM) 

P1020 


A(L1*ICOLUM)sO.O 

P1030 


DO 330 L*1*N 

PI 040 

330 

A(L1*L)*A(L1*L)-A(IC0LUM*L)*T 

P1050 


IF (M) 360*360*340 

P1060 

340 

DO 350 L*1*M 

P1070 

350 

B(L1»U=B(L1*L)-B(IC0LUM*L)*T 

P1080 

360 

CONTINUE 

P1090 

C 


PllOO 

c 

INTERCHANGE COLUMNS 

PlllO 

c 


P1120 


DO 390 1=1 *N 

P1130 


LsN*l-I 

P1140 


IF <INDEX(L*1)-INDEX(L*2) ) 370*390*370 

P1150 

370 

JR0W=1N0EX(L*1) 

P1160 


JC0LUM=IN0EX(L*2) 

P1170 


DO 380 Ksl*N 

P1180 


SWAP=A(K*JR0W) 

P1190 


A ( K ♦ JROW ) = A ( K * JCOLUM ) 

P1200 


79 



APPENDIX A - Continued 



A(K.JC0LUM)=SWAP 

P1210 

380 

Continue 

P1220 

390 

CONTINUE 

P1230 

400 

Return 

Pi 240 


End 

P1250- 
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APPENDIX A - Continued 


FUNCTION POLYEl (X,M,C) 0 10 

0 20 

This is a function used by subroutine lsqpol q so 

0 40 

Data big/ostttttttttt?/ q so 

Dimension c(m> o so 

IF CM-l) 30t40»10 0 70 

10 N=M-1 Q 80 

P0LYE1*C(1) Q 90 

Do 20 I=1»N Q 100 

20 POLYE1=X*POLYE1*C(I»1) Q llO 

RETURN 0 120 

30 P0LYE1=B16 Q 130 

RETURN Q 140 

40 POLYElsC(l) Q ISO 

RETURN Q 160 

End 0 170 
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APPENDIX A - Continued 


OVERLAY(RIN6,5.0) R 10 

Program critm r 20 

R 30 

This overlay computes critical MACH number R 40 

P 50 

Common /cord/ X<170),R(170)*C*RAVGfNN»N»OMIN»T(170).RB(170)*XBN R 60 

COMMON /EIVE/ CPO(170),MCR R 70 

Real mcr r bo 

Call mcrit <cpo»mcr) r 9o 

Return r loo 

End r 110- 
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APPENDIX A — Continued 


Subroutine mcrit <cpo»mcr) s lo 

S 20 

This subroutine computes critical MACH number of the nacelle BASEn S TO 

On the outer surface pressure S AO 

S 50 

Common /CORO/ x(170> »R(170) »C»RAVG*NN*N s 60 

Dimension cpouto) s ?o 

REAL MCR S 80 

AsO. S 90 

Do 30 JaltN S 100 

IF <CPO(J)) 10*30*30 S 110 

10 IF (A.EQ.O.) 60 TO 2o S 120 

IF (CPOU).LT.A) AsCPO(J) S 130 

Go TO 30 S 140 

20 AsCPO(J) S ISO 

30 CONTINUE S 160 

J=1 S 170 

McR»*999 S 1«0 

AO Cps(-2./MCR)*((l.-MCR)**1.5)»(S0RT(l.*MCR)> S 190 

IF (CP.LT.A) GO TO 50 S ?00 

MCRsMCR-.OOl 5 210 

JaJ*l S 220 

IF (J.6E.999) GO TO 50 S 230 

60 TO 40 S 240 

50 CONTINUE S 2 S 0 

RETURN s ?60 

End s 270- 
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APPENDIX A - Continued 


OvERLAY(RING,6*0) T 10 

Program masflow t 20 

C T 30 

c This overlay computes vi/vin by integrating the internal mass flow t ao 

C At some X STATION. THE X STATION IS PRESENTLY MIOCHORO BUT CAN T 41 

C BE CHANGED BY CHANGING CAROS T 170» T 320. T 650» T 970, T 42 

C T 50 

COMMON /CORO/ X(170> ,R(170) ,C,RAVG»NN.M.DMIN»T(170) ,RB(170) .XSI T 60 

COMMON /ONE/ ROUT(170) ♦RIN(170) *M.OMAX,XI (170) »XIO(170) *;AUPHO(l‘70) T 70 

1,ALPHOO(170»1) ,XBL(170) ,ALPT(170> ,ALPTH(170) »OROX(170| ♦OROxH(170) , T 80 

2LAM.IT.1TMAX T 90 

COMMON /TWO/ IP»CAM»6AM(170) »GQ(170) »ALPOO(170»1) T 100 

COMMON /SIX/ NRHO,RAD(170) »VIVIN,RXX»RBXI T 410 

DIMENSION RHOO(170) T 120 

Real lam.ksq.m t i30 

PRINT 150 T 140 

Print 13o t i5o 

PRINT 140 T 160 

XXX=C/2, T 170 

CALL MTLUP (XXX, RXX»2»NN, 170.1, IP, X,R) T 180 

Call mtlup (xxx,txx,2,nn,170,i,ip,x,T) t i9o 

RXXsRXX-(TXX/2.) T 200 

TeMP9*RXX/NRH0 T 210 

C T 220 

C COMPUTE ALL RAO(J) T 230 

C T 240 

00 20 Jsl.NRHO T 250 

IF (J.EQ.NRHO) 60 TO 10 T 260 

RAO(J)*RXX-( J*TEMP9) T 270 

GO TO 20 T 280 

10 RAD(J)=0, T 290 

20 CONTINUE T 300 

VIVIN=0, T 310 

XX=C/2, T 320 

call mtlup (XX,RBXI,2,NN,170,1,IP,X,RB) T 330 

IF (XX.LE.XBN) RBXisO. T 340 

TeSTsRBXI/RAVG T 350 

C T 360 

c Mass flow - do 9o t 370 

C T 380 

DO 90 J=l,NRHO T 390 

C T 400 

C COMPUTE THE RHOO(J) T 4l0 

C T 420 

IF (J.EQ.l) GO TO 30 T 430 

IF (J.EQ.NRHO) GO TO 40 T 440 

RH00(J)=(RA0(J-1)+RAD(J))/(2.*RAVG) T 450 

GO TO 50 T 460 

30 RhOO(J)=(RXX,RAD(J) )/(2.*RAVG) T 470 

GO TO 50 T 480 

40 RH00(J)=RA0(J-1)/(2,*RAVG) T490 

50 CONTINUE T 500 

IF (RHOO(J) .LE.TEST) GO TO 100 T 510 

C T 520 

c Velocities at rhoo(J> - do 6o t 530 

C T 540 

WXlJUOsO. T 550 

WRlJUOaO, T 560 

WX1QJU=0. T 570 

WR1QJU=0. T 580 
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WXOJU=0. T S90 

WROJUsO. T aOO 

WXBQJU=0, T MO 

WRBQJUsO, T 620 

XX»(C/N)*.25 T 1^30 

DO 60 I*1*N T 6^0 

DELXIsXKI) T 650 

KSQ*(^.*RH00(J) )/( (DELXI**2) ♦( (l.*RHOO( J) )«*2) ) T 660 

EKsEKA(KSQ) T (=<70 

SK=SKK(KSQ) T A80 

Call MTLUP (XX,RBXIt2,NN*170*l,IP,X.R8) T 690 

IF (XX.LE.XBN) RBXI=0. T 700 

TEMP1=(1./(2,*SQRT< (D£LXI**2) ♦< (1,*RH00(J) )**2) )))»(((( 2 . -KSQ )/( 1 . T 710 

1-KSQ) )*EK)-(2.*SK)-( (l./RHOO(J) )*(KSQ/(1.-KSQ) )*EK) ) T 7?0 

TeMP2=( ( (-1.)*DELXI)/(2.*RH00(J)*SQRT((0ELXI**2) ♦( (1,+RH00(J) )*»2) T 730 

1) ))*(<( (2. -KS0)/(1.-KS0) )»EK)-(2.*SK)> T 740 

tempt* ( 4,*0ELXI*EK)/( (S0RT( (DELXI**2)»( (1,+RH00( J) )**2) ) )<>( (OeLXI» T 760 

l*2)*(a.-RHOO(J))**2))) T 760 

TEMP8=(2,/(RHOO(J)*SQRT( (DELXI**2) ♦( (l.*RHOO(J) I **2 ) ) ) ) * ( ( ( 1 . - ( ( 2 . T 770 

1*RH00(J)*(RH00(J)-1.) )/( (DELXI**2) ♦( (l.-RHOO(J) )**2) ) ) )*EK)-SK) T 780 

TEMP11=RH00(J)/(2.*( (DElXI**2*RHOO(J)**2)**1.5) ) T 790 

TEMP12s(DELXn/(2.*( (DELXI**2*RH00( J)**2)**l .5) ) T «00 

WX1JU0*WXIJU0*TEMP1*ALPH0{I)*GAM<I)*<1./(2.*3.14159) ) T MO 

Wr1JU0sWR1JU0+TEMP2*ALPH0(I)*GAM(I)*(1./<2.*3.14159) ) T 820 

WX10JUsWX10JU*ITENP1*G0(I)*«XI(1)-XH2))/{2.*3.1A159)) T 830 

WR1QJUsWR1QJU*(TEMP2»G0(1)*(XI{1)-XI(2))/(2.*3.14159) ) T «40 

WXOJUsWXQJU-(TEMp-/»ALPT(I)*(XI (l)-XI (2) )/<2.*3.14159) ) T 860 

WRQJU=WR0JU+(TEMP8*ALPT(I)*(XI(1)-XI(2))/(2.*3.14169) ) T B60 

WXB0JU=WXB0JU+TEMP12*(RBXI/RAVG)*0R0X(I)*(XI(1)-XI(2) ) T 870 

WRBQJU=WRB0JU+TEMP11*(RBXI/RAVG)*0RDX(I)«(XK1)-XI (2) ) T 880 

XX*XX*(C/N) T 890 

60 CONTINUE T oOO 

WXlJUOIsWXlJUO T qio 

WxiQJUIaWXlOJU T 920 

VTI=U.*(WX1JU0I+WX10JUI*WXQJU+WXBQJU)/(1.-M**2) )**2 T 930 

C T 940 

C SUMMATION' FOR MASS FLOW - ADO ONE TERM EACH TIME THRU 00 9o T 950 

C T 960 

XX=C/2. T 970 

Call mtlup (xx,rbxi,2.nn.ito»i,ip,x,rb) t 980 

IF (XX.LE.XBN) RBXI=0. T 990 

IF (J.EQ.l) GO TO 70 TlOOO 

TEMP10=S0RT(VTI)*( (RAD( J-1)*»2-RA0(J)»*2)/(RXX**2-RBXI**2) ) TlOlO 

60 TO 80 T1020 

70 TeMP10sSQRT(VTI)*( (RXX**2-RAD(J)**2)/(RXX**2-WBXI*<*2) ) T1030 

80 CONTINUE Tlo^O 

V1VIN=VIVIN+TEMP10 TloSO 

WRITE (6»120) J»RAD< J) «RH00( J) .WXlJUO»WRlJu6»WXlCIJU»WRlOJU,WXfoJU,W T1060 

IRQJU T1070 

write (6.110) VTI.WIVIN.TEMPIO.WXBOJU.WRBOJU Tl0«0 

90 CONTINUE T1090 

100 CONTINUE TllOO 

RETURN Til 10 

C T1120 

C T1130 

110 Format (13x.s(2x,eii.4)/) tii40 

120 format (1X»I3.3(2X,F10.5) »6(2X,EU.4) ) T1150 

130 FORMAT ( IX. 1HJ.9X. 3HRAn»8X»4HRHP0.6X»6HWXl JUO . 7X . 6HWR 1 JUO * 7X . 6HWX 1 T1160 

IQJU.7X.6HWRI0JU.8X.5HWXQJU.8X.SHWRQJU) Til 70 

140 FORMAT (23X .3HVTI . 8X.5HVIV1N. 7X .6HTEMP1 0.7X . 6HWXB0JU. 7X . 6MwRBU JU/ ) T1180 
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APPENDIX A - Concluded 


ISO FORMAT (lHlf21HMASS FLOW COMPUTATION* 10X»49HALL PERTURB VEl AND RA T1190 
ID ARE FOR EQUIV INCOMP 6E0M//> T1200 

End T1210- 
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APPENDIX B 


DERIVATION OF EQUATIONS (13) AND (23) 


The velocity induced in the x-direction at ?q>Pq by a distribution of vortex rings 
on a cylinder of radius p = 1 and of chord length c is found by use of equation (21). 
The result is 


J-l/\ 4 


2 - 




a - k2 pQ 1 



E - 2H 


d| (Bl) 


Replacing the modulus k by its definition in equation (11) and rearranging the result, 
equation (Bl) becomes 


!hdVo),j. rVx 

27rJ-1/X 


E - K 




■ 1 /x «„(«)ro® 2(p - 1) 




E 




d^ 


(B2) 


Equation (B2) must now be evaluated at Pq = 1 to obtain equation (13). Considering the 
first integral, at Pq = 1 and | = ^q, a singularity exists in the term containing the ellip- 
tic integral of the first kind K. To investigate this singularity, examine the integral in a 
small region extending e on each side of 


lim 

e-0 


ji. 

2tt j t _ 6 
^o 


tto(i) yo(i)K d{ 


~^o(^o)^o(^o) 

27T 


lim 

e-0 



Kd^ 


(B3) 


where the terms multiplying K are approximated by setting | The elliptic inte- 

gral K can be expanded in a series as a function of the complimentary modulus k’ 

(ref. 23) in the following form 



((k')^ < l) (B4) 
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APPENDIX B - Continued 


where 



and Pochhammer's symbol is defined by 



and is the digamma function (ref, 24). By using equation (B4) and letting e — 0, and 
therefore k' -*• 0, it can be shown that the integral in equation (B3) goes to 0. Thus the 
principal value of the first integral in equation (B2) is found by treating it as an ordinary 
integral. 

For Pq = 1, the value of the second integral in equation (B2) can be determined as 
follows. For ^ the integrand is 0 when = 1. Hence the limit must be found as 
^ simultaneously with Pq — 1 



where E “ 1 because k » 1, and again the first term in the integrand is approximated 
by setting ^ By defining a new variable 



and requiring p^ -► 1 faster than e -• 0 for a nontrivial solution, it is found that 
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APPENDIX B - Concluded 



so that 


to ao({)ro({) 2 (Pq-1)e 

+ (p„ - 1) 


,B5) 

2 2 A| 


where = ^o(^o)/^^* substituting equation (B5) into equation (B2) and setting 


Uy(^o>^) ^ %(4)yo(g) 


E - K ^ %(^o) ^o(^o) 

. t \2 . .1 2 


m-m 


Writing equation (B6) in summation form yields equation (13). Equation (23) is derived by 
a procedure analogous to the preceding derivation. 
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